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LIGA  is  a  relatively  new  technology  used  to  pattern  high  aspect  ratio  microelectromechanical  systems  (HARMEMS)  in  a 
resist  material,  using  X-ray  radiation.  Resist  materials  used  in  LIGA  typically  have  conduction  and  thermal  expansion 
properties  that  are  very  different  from  those  of  the  substrate  that  supports  them  during  the  exposure  process,  and  thermal 
deformations  may  limit  the  ability  of  the  LIGA  process  to  reproduce  patterns  accurately  in  the  resist.  A  knowledge  of  the 
temperature  distributions  in  the  resist  and  substrate  will  facilitate  the  study  of  thermal  deformations  and  their  effects  on  the 
manufacturing  process. 

The  solution  of  analytical  models  of  the  temperature  distributions  in  the  resist-substrate  system  are  presented.  The  primary 
models  presented  are  a  one-layer,  two-dimensional  model,  and  a  two-layer  model  in  one  dimension.  Boundary  conditions  are 
developed  based  on  current  practices  used  in  the  LIGA  process. 

Subjects  relating  to  the  evaluation  of  the  solutions  are  discussed,  including  the  characteristics  of  series  solutions  and  the 
development  of  computer  programs  to  handle  calculations. 

The  results  of  some  simple  temperature  measurement  experiments  performed  at  the  Center  for  Advanced  Microstructures 
and  Devices  (CAMD)  are  presented,  along  with  a  discussion  of  the  relative  merits  of  experimental,  computational  and 
analytical  methods  of  analysis. 
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ABSTRACT 


LIGA  is  a  relatively  new  technology  used  to  pattern  high  aspect  ratio  microelectromechanical  systems 
(HARMEMS)  in  a  resist  material,  using  X-ray  radiation.  Resist  materials  used  in  LIGA  typically  have 
conduction  and  thermal  expansion  properties  that  are  very  different  from  those  of  the  substrate  that  supports 
them  during  the  exposure  process,  and  thermal  deformations  may  limit  the  ability  of  the  LIGA  process  to 
reproduce  patterns  accurately  in  the  resist.  A  knowledge  of  the  temperature  distributions  in  the  resist  and 
substrate  will  facilitate  the  study  of  thermal  deformations  and  their  effects  on  the  manufacturing  process. 

This  thesis  presents  the  solution  of  analytical  models  of  the  temperature  distributions  in  the  resist- 
substrate  system.  The  primary  models  presented  are  a  one-layer,  two-dimensional  model,  and  a  two-layer 
model  in  one  dimension.  Boundary  conditions  are  developed  based  on  current  practices  used  in  the  LIGA 
process. 

Subjects  relating  to  evaluation  of  the  solutions  are  discussed,  including  the  characteristics  of  series 
solutions  and  the  development  of  computer  programs  to  handle  the  calculations. 

The  results  of  some  simple  temperature  measurement  experiments  performed  at  the  Center  for 
Advanced  Microstructures  and  Devices  (CAMD)  are  presented,  along  with  a  discussion  of  the  relative 
merits  of  experimental,  computational,  and  analytical  methods  of  analysis. 


APPROVAL  FOR  SCHOLARLY  DISSEMINATION 


The  author  grants  to  the  Prescott  Memorial  Library  of  Louisiana  Tech  Univeisity 
the  right  to  reproduce,  by  ^propriate  methods,  upon  request,  any  or  all  portions  of  this 
Thesis.  It  is  understood  that  "proper  request"  consists  of  the  agreement,  on  the  part  of  the 
requesting  party,  that  said  reproduction  is  for  his  personal  use  and  that  subsequent 
reproduction  will  not  occur  without  written  ^proval  of  the  author  of  this  Thesis.  Further, 
any  portions  of  the  Thesis  used  in  books,  p^ers,  and  other  works  must  be  appropriately 
referenced  to  this  Thesis. 

Finally,  the  author  of  this  Thesis  reserves  the  right  to  publish  freely,  in  the 
literature,  at  any  time,  any  or  all  portions  of  this  Thesis. 


Author _ -4  G 

Date _ *?/io  /q(s 


GS  Form  14 

(2/88) 


TABLE  OF  CONTENTS 


ABSTRACT  .  iii 

LIST  OF  FIGURES  .  vii 

LIST  OF  TABLES  .  viii 

LIST  OF  SYMBOLS  .  ix 

Chapter 

1.  INTRODUCTION  .  1 

2.  RELATED  RESEARCH  .  4 

3.  PROBLEM  DEHNITION  .  8 

3.1.  Physical  System  .  8 

3.2.  Two-Dimensional  Model  .  12 

3.3.  Two-Layer  Model  .  16 

3.4.  Problems  With  the  Two-Layer,  Two-Dimensional  Model  .  18 

4.  MATHEMATICAL  SOLUTION  .  21 

4.1.  One-Layer  Solution  .  21 

4.2.  Steady-State  Solution  .  31 

4.3.  Two-Layer,  Time-Dependent  Solution  .  35 

5.  ANALYSIS  OF  SOLUTION  .  47 

5.1.  Series  Solutions  .  47 

5.2.  Computer  Programs  .  50 

5.3.  Test  Cases  .  54 

5.4.  Experimental  Results  .  58 

5.5.  Conclusions  .  62 


V 


6.  CONCLUSIONS  AND  RECOMMENDATIONS  .  64 

6. 1  Analytical  Models  .  64 

6.2  Experimental  Analysis  .  65 

6.3  Recommendations  .  66 

APPENDDC:  COMPUTER  CODE  .  67 

BIBLIOGRAPHY  .  87 

VITA  .  89 


VI 


LIST  OF  HGURES 


Page 

1.1.  Typical  LIGA  exposure  station  design  .  3 

3.1.  Physical  model  for  the  one-layer  problem  .  1 3 

3.2.  Physical  model  for  the  two-layer  problem  .  1 7 

4. 1 .  Temperature  rise  for  generation  part  of  one-layer  solution  .  28 

4.2.  Temperature  for  multiple  passes  of  the  X-ray  beam  at  center  of  resist  .  31 

4.3.  Temperature  distribution  for  steady  state  test  of  two-layer  solution  .  44 

5.1.  Convergence  of  one-layer  solution  in  x  direction  at  y  =  L  .  48 

5.2.  Convergence  of  one-layer  solution  in  y  direction  at  the  center  of  the  resist  .  49 

5.3.  Convergence  of  one  and  two-layer  solutions  .  49 

5.4.  Graph  of  the  eigenvalue  function  for  one-layer  problem  .  51 

5.5.  Determinant  of  the  eigenvalue  matrix  for  two-layer  problem  .  52 

5.6.  Temperatures  from  test  of  the  substrate  temperature  gradient  .  56 

5.7.  Temperatures  with  respect  to  time  for  one-layer  geometry  .  58 

5.8.  Wafer  used  in  experimental  test  of  heat  transfer  .  59 

5.9.  Temperature  ranges  for  the  first  exposure  test  .  61 

5.10.  Temperature  ranges  for  the  second  exposure  test  .  61 

5.11.  Temperature  ranges  for  the  third  exposure  test  .  62 


LIST  OF  TABLES 

Table  Page 

3. 1 .  Dosage  profiles  for  a  typical  wafer  at  CAMD  .  1 5 

4.1.  Constants  used  for  generation  solution  .  28 

4.2.  Constants  for  test  of  the  steady-state  solution  .  34 

4.3.  Boundary  point  temperatures  for  steady-state  solution  . 35 

4.4.  Constants  for  test  of  two-layer  solution  with  no  generation  .  42 

4.5.  Boundary  point  temperatures  for  steady-state  test  .  43 

4.6.  Constants  for  test  of  the  resist  generation  term  .  45 

4.7.  Temperatures  (in  °C)  from  tests  of  the  two-layer  solution  .  45 

5. 1 .  Constants  for  test  of  the  substrate  temperature  gradient  .  56 

5.2.  Constants  for  calculation  of  multiple  pass  temperatures  .  57 


LIST  OF  SYMBOLS 


a  Width  of  wafer  in  x  direction 

A  Area  Coefficient  for  one  dimensional  test  solution,  area 

Af  fjj  Eigenfunction  coefficient 

B  Coefficient  for  one-dimensional  test  solution 

^  Eigenfunction  coefficient 

c  Specific  heat 

C  Coefficient  for  one-dimensional  test  solution 

C, .  Coefficient  for  two-layer,  steady-state  solution 

d  Distance  between  mask  and  resist 

D  Coefficient  for  one-dimensional  test  solution 

D.  Coefficient  for  two- layer,  steady-state  solution 

F(x,y)  Initial  condition 

g  Generation  function 

h\  Heat  transfer  coefficient  at  the  resist  surface 

/t2  Contact  resistance  coefficient 

Heat  transfer  coefficient  at  the  substrate  surface 
lij  divided  by  the  appropriate  -  a  simplification 
k  Thermal  conductivity 

L|  Thickness  of  the  resist 

L2  Thickness  of  wafer  (resist  and  substrate) 

N  Orthogonality  constant 

p  Pressure 

q  Heat  transfer 

r  General  spatial  coordinate 

t  Time 

T  Temperature  rise 

7]  Ambient  temperature  at  the  substrate  surface 

Ty  Ambient  temperature  at  the  resist  surface 

u  Gas  velocity 


IX 


V 

Velocity  of  the  X-ray  beam  with  respect  to  the  resist 

Wafer  velocity  for  gas  velocity  solution 

w 

Width  of  the  X-ray  beam  in  the  x  direction 

W 

Irradiance 

X 

Spatial  variable  along  the  length  of  the  wafer 

X{x) 

X-  dependent  part  of  separated  solution 

y 

Spatial  variable  through  the  depth  of  the  wafer 

ny) 

j-dependent  part  of  separated  solution 

Greek  Alohabet 

a 

Thermal  diffusivity 

p 

X  direction  eigenvalues 

Y 

Simplification  variable;  defined  in  Eq.  (4.69) 

Eigenvalues  for  gas  velocity  solution 

n 

Simplification  variable;  defined  in  Eq.  (4.69) 

e(0 

Time-dependent  part  of  separated  solution 

A 

y-dircction  eigenvalues 

P 

Absorption  coefficient  or  viscosity 

V 

Kinematic  viscosity 

P 

Density 

H' 

Steady-state  part  of  two-dimensional  solution 

Subscripts 

I 

Referring  to  the  resist 

2 

Referring  to  the  substrate 

i 

Layer  1  or  2  (resist  and  substrate,  respectively) 

1 

Lower  surface  (bottom  of  resist  or  substrate) 

m,  n 

One  of  an  infinite  series,  y  or  jc  direction,  respectively 

u 

Upper  surface  (top  of  resist) 

d 

Decay  solution 

CHAPTER  I 


INTRODUCTION 

The  LIGA  process  (a  German  acronym  for  Lithographic,  Galvanoformung,  Abformung)  is  a  relatively 
new  technology  used  to  pattern  high  aspect  ratio  microelectromechanical  systems  (HARMEMS)  in  a  resist 
material.  The  process  is  similar  to  optical  lithography  used  for  the  manufacture  of  microchips.  X-ray 
radiation  passes  through  a  mask  to  a  resist  material,  which  is  chemically  altered  by  the  incident  radiation. 
The  mask  is  patterned  with  an  absorbing  material,  and  a  negative  of  this  pattern  is  transferred  to  the  resist. 
A  variety  of  substances  have  been  used  for  resist  materials,  but  the  most  common  is 
Polymethylmethacrylate,  or  PMMA.  The  resist  must  be  supported  by  a  substrate;  the  choice  of  materials 
for  this  component  is  more  flexible,  and  depends  on  the  application.  The  bonded  resist  and  substrate  will 
be  referred  to  as  a  wafer.  In  many  cases  a  layer  of  some  other  material,  such  as  a  metal  base,  is  applied  to 
the  substrate  before  addition  of  the  PMMA,  to  facilitate  later  processing.  A  metal  layer  allows 
electroplating  of  metal  into  the  gaps  that  remain  after  the  resist  has  been  chemically  treated  to  remove  the 
material  which  was  altered  by  the  X-ray  radiation. 

The  advantage  LIGA  has  over  conventional  lithography  is  the  ability  to  produce  structures  with  far 
higher  aspect  ratios.  While  the  lateral  dimensions  of  the  structures  may  be  accurate  to  a  fraction  of  a 
micron,  their  height  can  be  hundreds  of  microns.  This  flexibility  has  allowed  the  design  of  a  variety  of 
functional  devices  and  systems,  including  such  things  as  sensors,  microactuators,  various  fluidic 
components,  and  microoptics. 

A  diagram  of  a  typical  LIGA  exposure  station  is  shown  in  Fig.  1.1.  During  the  exposure  process,  the 
X-ray  beam  is  held  steady  while  the  carrier,  which  holds  the  mask  and  wafer,  is  moved  up  and  down.  The 
beam  is  typically  rectangular  in  cross  section,  with  a  height  on  the  order  of  a  centimeter,  and  several 
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centimeters  wide.  Sweeping  the  wafer  through  the  beam  is  done  to  give  an  even  exposure  of  the  resist  and 
to  avoid  damaging  one  part  of  the  resist  by  overexposure. 

Synchrotrons  are  used  to  produce  the  X-rays  for  the  LIGA  process.  Sources  such  as  the  one  at  the 
Louisiana  State  University  Center  for  Advanced  Microstructures  and  Devices  (CAMD)  in  Baton  Rouge  can 
produce  high  flux,  long  wavelength  radiation,  which  requires  exposure  times  of  only  a  few  minutes.  This  is 
advantageous  because  of  the  high  expense  of  using  these  facilities.  However,  the  heat  generated  within  the 
resist  and  substrate  by  the  incident  radiation  may  cause  problems. 

Thermal  stresses  are  generated  within  the  resist  material  itself,  because  only  part  of  the  resist  is  being 
exposed  to  radiation.  Also,  resist  materials  typically  have  very  different  conduction  and  thermal  expansion 
properties  than  the  substrate.  Thermal  stresses  within  the  resist  and  substrate  may  limit  the  ability  of  the 
LIGA  process  to  accurately  reproduce  patterns  in  the  resist.  Thermal  deformations  may  cause  the  resist  to 
move  with  respect  to  the  mask,  causing  incomplete  exposure.  In  addition,  the  resist  may  move  during 
cooling  after  the  exposure  process  is  finished,  shifting  away  from  the  pattern  that  was  imprinted  on  the 
heated  structure.  Thermal  stresses  may  also  contribute  to  bonding  problems  between  the  resist  and 
substrate;  the  two  materials  often  separate  during  the  exposure  or  development  processes. 

Knowledge  of  the  temperature  distributions  in  the  resist  and  substrate  will  facilitate  the  study  of 
thermal  deformations,  and  hopefully  allow  researchers  to  minimize  the  effects  of  those  deformations  on  the 
manufacturing  process.  This  thesis  presents  the  solution  of  analytical  models  of  the  temperature 
distributions  in  the  resist-substrate  system.  The  primary  models  presented  are  a  one-layer,  two-dimensional 
model,  and  one  with  two  layers  and  one  dimension.  Boundary  conditions  are  developed  based  on  current 
practices  used  in  the  LIGA  process,  while  allowing  for  flexibility  in  analyzing  the  effects  of  different 
configurations  and  exposure  conditions. 


3 


Fig  1.1.  Typical  LIGA  exposure  station  design 


The  original  goal  of  this  thesis  was  to  model  the  resist-substrate  system  with  a  two-layer,  two- 
dimensional  solution  of  the  heat  transfer  equations.  No  analytical  method  was  found  to  accomplish  this.  A 
number  of  authors  have  indicated  that  this  type  of  solution  is  possible,  notably  Osisik  (1980).  A  solution 
was  derived  using  the  method  outlined  by  Ozisik  (1980),  and  was  found  to  be  incorrect;  a  lengthy 
investigation  revealed  that  the  method  itself  is  at  fault.  The  reasons  for  this  are  outlined  in  Chapter  3. 


CHAPTER  2 


RELATED  RESEARCH 

High  flux  X-ray  lithography  is  a  relatively  new  technology  which  can  be  praeticed  at  only  a  few 
facilities  in  the  world.  So,  despite  the  fact  that  excessive  heat  generation  and  heat  induced  deformations  are 
common  problems  when  very  high  fluxes  are  used,  little  has  been  published  on  these  subjects. 

Some  preliminary  research  into  resist  heating  was  done  by  Ameel  et  al.  (1994).  Their  work  presents 
three  simple  analytical  models  of  resist  temperature  distributions.  The  models  are  based  in  a  cylindrical 
coordinate  system,  and  consider  only  the  resist  material;  the  difference  between  them  is  the  boundary 
condition  at  the  bottom  of  the  resist,  which  approximates  the  interface  with  the  substrate.  An  attempt  was 
made  to  find  minimum  and  maximum  temperature  rises,  establishing  a  range  of  temperatures  one  might 
expect  to  see  in  practice.  For  the  specific  case  studied,  this  range  was  from  2.3“C  to  1 1.8°C.  The  large 
range,  along  with  the  simplicity  of  the  model  and  boundary  conditions,  limits  the  utility  of  these  models. 

A  large  amount  of  the  X-ray  radiation  used  in  LIGA  is  absorbed  by  the  mask.  It  is  much  thinner  and 
has  less  volume  than  the  resist-substrate  system,  and  is  harder  to  cool,  so  temperature  rises  in  the  mask  are 
generally  expected  to  be  higher.  For  this  reason,  some  effort  has  been  put  into  determining  heat  generation 
and  deformation  in  LIGA  masks,  while  less  of  this  type  of  work  has  been  done  on  the  resist  and  substrate. 

A  finite  element  analysis  of  temperature  distributions  in  X-ray  irradiated  masks  for  the  LIGA  process 
has  been  reported  by  Feiertag  et  al.  (1994).  Actual  surface  temperature  measurements  were  used  to  confirm 
the  predictions  of  the  model.  The  mask,  like  the  resist,  is  composed  of  two  layers  of  material  with  differing 
thermal  properties.  However,  the  two  situations  are  fundamentally  different;  thick  mask  carriers  increase 
the  necessary  exposure  times,  limiting  their  ability  to  act  as  heat  sinks;  thick  resist  substrates  may  be  used  to 
absorb  heat  and  reduce  deformation  of  the  resist-substrate  structure. 
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This  work  by  Feiertag  et  al.  (1994)  is  useful  for  its  analysis  of  the  conditions  in  the  gap  between  the 
mask  and  the  resist;  this  is  an  important  aspect  of  the  resist  model.  The  authors  assumed  negligible  heat 
loss  due  to  convection  and  radiation  from  the  surface  of  the  mask  and  carrier.  Their  arguments  are  taken 
into  consideration  in  the  determination  of  appropriate  boundary  conditions  for  the  present  model. 

Vladimirsky  et  al.  (1989)  present  an  analysis  of  heating  and  thermal  deformations  in  an  X-ray 
irradiated  mask,  with  both  a  theoretical  model  and  experimental  surface  temperature  measurements.  The 
experimental  results  did  not  compare  well  with  the  calculations  for  most  of  the  cases  considered,  but  the 
authors  made  an  attempt  to  explain  these  discrepancies.  Variation  of  the  mask  temperature  with  the  spacing 
of  the  mask-resist  gap  was  documented.  The  temperature  data  and  observations  about  heat  transfer  in  the 
gap  have  applications  to  the  determination  of  appropriate  boundary  conditions  for  the  project  reported 
herein. 

The  application  of  mathematical  solutions  to  the  solution  of  conduction  problems  in  multiple 
dimensions  was  begun  in  the  1930s.  Only  the  simplest  solutions  may  be  calculated  without  the  use  of  a 
computer,  and  the  subject  generated  little  interest  until  the  1960s,  judging  by  the  amount  of  work  published. 
Osisik  (1980)  published  the  first  edition  of  his  Heat  Conduction  text  in  the  mid-60s;  it  is  the  most  recent 
edition  which  is  used  as  a  reference  for  this  work,  for  almost  all  of  the  analytical  conduction  solutions. 
Many  of  the  papers  published  since  the  first  edition  of  Ozisik  (1964)  list  it  as  a  reference.  The  general 
solution  in  Ozisik  (1980)  provided  the  framework  on  which  the  two-dimensional,  two-layer  solution  was 
based.  The  problems  with  this  general  solution  are  presented  in  Chapter  3.  Cobble  (1970)  also  presents  a 
detailed  analysis  of  the  problem  of  multiple-layer  conduction,  though  his  work  is  limited  to  one  dimension. 

The  problem  of  heat  generation  in  a  multilayered  system  has  a  number  of  engineering  applications, 
and  many  papers  have  been  published  on  this  topic.  One  situation  closely  related  to  the  present  problem  is 
that  of  laser  heating  of  thin  films,  a  common  problem  in  laser  optical  systems.  The  layers  represent  optical 
coatings,  which  rest  on  a  mirror,  modeled  as  an  infinite  solid.  These  studies  generally  use  cylindrical 
coordinates,  but  have  some  similarity  to  the  problem  being  studied  here.  Many  simplifications  and 
assumptions  are  made  in  these  studies  that  would  not  be  applicable  to  the  LIGA  system.  The  authors  are 
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generally  more  concerned  with  reflection  and  absorption  in  the  film  layers  than  with  the  conduction  of  heat 
between  them. 

Cole  and  McGahan  (1993)  obtained  a  partly  analytic  solution  for  laser  heating  of  a  multilayer  system. 
They  used  local  Green's  functions  for  each  layer,  along  with  Fourier  and  Hankel  integral  transforms,  but  the 
method  required  numerical  integration  to  get  a  temperature  distribution.  The  model  took  into  account  the 
contact  resistance  between  layers,  and  heat  generation  within  each  layer.  Though  the  author’s  model  was 
based  on  a  cylindrical  coordinate  system,  the  procedure  could  probably  be  applied  to  the  LIGA  problem; 
however,  the  numerical  part  is  unavoidable,  and  not  in  accordance  with  the  object  of  this  thesis. 

The  problem  of  laser  heating  in  multiple  layers  was  also  studied  by  Kant  (1988).  Laplace  and  Hankel 
transforms  were  used  to  construct  a  two-dimensional  model  in  the  cylindrical  coordinate  system,  and  the 
author  indicates  that  the  procedure  is  applicable  to  any  number  of  layers.  As  in  Cole  and  McGahan  (1993), 
numerical  methods  are  needed  to  find  the  inversions  of  the  transforms.  The  model  is  very  similar  to  the  one 
used  in  Cole  and  McGahan  (1993),  and  in  work  by  Madison  and  McDaniel  (1989).  The  latter  authors 
studied  the  same  basic  system,  though  they  did  use  an  exponential  generation  term  similar  to  the  one  used  in 
this  thesis. 

El-Adawi  et  al.  (1995)  studied  laser  heating  in  a  two-layer  system.  Their  model  was  one-dimensional, 
and  the  only  significant  difference  between  it  and  the  two-layer  solution  in  the  present  work  is  that  the  heat 
generation  was  introduced  through  a  constant  flux  boundary  condition,  as  opposed  to  volumetric  heating. 
The  solution  was  obtained  with  Laplace  transforms;  this  method  is  very  different  form  the  one  used  here, 
but  no  less  complicated.  The  authors  made  no  mention  of  the  extension  of  the  problem  to  more  than  one 
dimension. 

Mikhailov  (1973)  introduced  a  finite  integral  transform  and  inversion  for  the  solution  of  the  diffusion 
equation  in  an  arbitrary  region  with  coupled  boundary  conditions.  The  derivations  were  detailed  but  very 
general.  The  solution  involves  the  derivation  of  eigenfunctions  and  eigenvalues  from  a  corresponding 
Sturm-Louiville  problem  in  two  regions.  Mikhailov  notes  that  although  the  eigenfunctions  will  be  different 
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for  the  two  regions,  they  must  share  the  same  eigenvalues.  This  requirement  is  related  to  the  problems 
encountered  in  finding  a  two  layer  solution  for  more  than  one  dimension  (see  Chapter  3). 

Sareen  and  Gidaspow  (1974)  studied  mass  diffusion  in  a  two-layer  system.  The  mathematical  model 
was  almost  identical  to  the  model  proposed  for  this  thesis,  but  with  one  major  assumption  -  that  the 
diffusivity  is  the  same  for  both  layers  in  the  lengthwise  direction.  For  the  LIGA  model,  this  would  mean 
assuming  that  the  thermal  diffusivity  is  the  same  for  the  resist  and  substrate  (at  least  in  the  direction 
tangential  to  the  surface),  when  in  fact  they  differ  by  several  orders  of  magnitude.  This  assumption  and  its 
implications  are  discussed  in  more  detail  in  Chapter  3. 

The  analytical  work  presented  here  is  basically  a  continuation  of  work  done  by  Ameel  et  al.  (1994). 
That  work  was  the  only  one  found  which  directly  addressed  the  problem  of  X-ray  heating  of  a  LIGA  resist. 
The  multiple-layer  problems  described  above  are  mathematically  similar  to  some  of  the  problems  developed 
herein,  though  the  physical  situation  is  very  different.  The  works  dealing  with  laser  heating  tended  to 
concentrate  on  the  effects  of  reflection  and  diffraction,  and  mixed  analytical  and  numerical  solution 
procedures. 

The  two-dimensional,  two-layer  solution  which  was  to  be  the  object  of  the  present  work  was  not  found 
in  the  literature  in  any  form;  a  number  of  authors  mentioned  that  such  a  solution  is  possible,  giving  Ozisik 
(1980)  as  a  reference,  but  nobody  (including  Ozisik)  developed  a  specific  solution  to  that  problem. 


CHAPTER  3 


PROBLEM  DEFINITION 

The  models  used  to  predict  the  temperature  field  should  be  based  as  closely  as  possible  on  the  physical 
system,  consisting  of  the  bonded  resist  and  substrate  (the  wafer),  the  carrier  which  holds  the  mask  and 
wafer,  and  the  physical  surroundings.  It  is  also  important  to  determine  appropriate  values  for  heat  transfer 
coefficients  and  the  thermodynamic  properties  such  as  thermal  conductance.  Many  of  these  values  are 
easily  found,  but  others  require  a  number  of  assumptions  and  calculations.  In  this  chapter,  the  physical 
situation  will  be  examined,  and  appropriate  models,  boundary  conditions,  and  constants  determined. 
Although  the  originally  envisioned  problem  was  not  able  to  be  solved  in  its  complete  form,  some  simpler 
models  will  be  developed  which  may  give  useful  data. 


3.1.  Physical  System 

The  wafer  is  flat  and  may  be  either  circular  or  rectangular.  The  heat  transfer  within  the  wafer  may  be 
easily  represented  by  general  heat  conduction  equations; 


g(r) 


dT 


a  =  ■ 


(3.1) 


k  (X  dt  pc 

where  T  is  temperature,  g  is  a  generation  function,  t  is  time,  k  is  thermal  conductivity,  a  is  thermal 
diffusivity,  and  r  is  a  general  spatial  variable.  The  form  of  the  Eq.  3.1  assumes  thermal  conductivity  to  be 
constant;  this  is  a  good  assumption,  given  the  small  temperature  ranges  which  will  be  shown  to  occur. 


One  of  the  boundary  conditions  on  Eq.  3.1  depends  on  the  carrier,  which  supports  the  wafer  during  the 
exposure  process,  and  its  relation  to  the  substrate.  Because  stability  is  important  for  maintaining  proper 
alignment  of  mask  and  wafer,  the  carrier  is  very  massive  in  comparison  the  resist.  The  wafer  is  clamped 
directly  on  the  metal  surface  of  the  carrier,  substrate  side  down,  with  the  mask  assembly  in  front  of  it. 
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There  will  be  some  resistance  to  conduction  at  the  interface  between  substrate  and  carrier,  with  the  carrier 
acting  essentially  as  an  infinite  heat  sink  at  constant  temperature. 

The  boundaiy  condition  on  Eq,  3.1  that  is  opposite  the  carrier  side  of  the  wafer  depends  on  conditions 
in  the  exposure  station.  Jara-Almonte  (1995)  has  indicated  that  chamber  conditions  vary  with  the  situation 
and  researcher,  but  that  he  uses  a  mixture  consisting  of  about  98%  helium  at  a  pressure  of  approximately  10 
kPa,  or  one  tenth  of  an  atmosphere,  with  the  gas  at  approximately  room  temperature.  The  carrier  assembly, 
with  the  wafer  inside,  is  moved  up  and  down  through  this  gas.  This  surface  may  be  water  cooled,  but  the 
same  boundary  condition  can  be  used  for  that  situation,  with  a  change  in  constants.  The  resist  surface, 
however,  is  separated  from  the  mask  by  a  space  typically  on  the  order  of  tens  of  micrometers  (microns). 
This  gap  will  allow  for  very  little  gas  flow  through  the  space,  and  the  oscillation  of  the  system  will  act  to 
further  reduce  constant  circulation. 

3.1.1.  Micro  Effects 

If  the  scale  of  a  system  is  sufficiently  small,  the  relations  commonly  used  for  determination  of  physical 
phenomenon  may  be  inaccurate.  This  is  because  the  traditional  relations  treat  matter  as  a  continuous 
medium.  Of  course,  this  is  a  simplification;  matter  is  made  up  of  atoms,  and  when  dealing  with  very  small 
geometries  interactions  between  the  individual  atoms  may  become  important.  A  high  vacuum,  in  which  the 
molecules  are  widely  spaced,  may  also  necessitate  consideration  of  individual  molecules.  Because  the  gap 
between  the  resist  and  mask  has  both  of  these  characteristics,  a  detailed  analysis  of  this  region  is  needed. 

First,  consider  an  estimate  of  the  mean  free  path  of  helium  atoms  under  the  conditions  given  in  the 
previous  section.  This  quantity  is  easily  calculated,  given  the  right  properties;  Keesom  (1942)  gives  some 
formulas  and  representative  values.  At  one  tenth  of  an  atmosphere  pressure  and  300  K,  the  mean  free  path 
is  about  2.0  x  10  ^mm  .  The  distance  between  the  mask  and  the  resist  should  be  taken  as  the  characteristic 
dimension.  This  distance  is  generally  less  than  0.5  mm  (Feiertag,  1994),  and  the  temperature  of  the  gas  in 
the  gap  will  increase  slightly  during  exposure.  However,  the  ratio  of  mean  free  path  to  the  characteristic 
dimension  (the  Knudsen  number)  should  still  be  on  the  order  of  about  0.01.  This  is  not  large  enough  to 
cause  a  significant  deviation  from  the  classical  material  properties  used  in  following  discussions. 
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The  resist  layer,  though  thin  by  most  standards  (typically  less  than  half  a  millimeter),  is  very  thick 
when  compared  to  the  spacing  of  the  molecules  of  which  it  is  made,  and  will  be  governed  by  macro  scale 
relations.  However,  there  are  often  one  or  more  layers  between  the  mask  and  the  resist,  some  of  which  may 
be  thin  enough  to  have  nonclassical  conduction  characteristics.  For  instance,  a  layer  of  metal  is  often  plated 
onto  substrate  materials  before  application  of  the  resist,  to  act  as  a  base  for  later  electroplating  processes, 
and  antireflective  coatings  may  be  placed  on  the  resist  to  keep  reflected  radiation  from  affecting  the  resist. 
Though  the  method  used  for  the  two-layer  model  here  is  applicable  to  any  number  of  layers,  additional 
layers  would  greatly  complicate  it.  Little  information  concerning  the  details  of  these  coatings  and  exposures 
has  been  published;  it  would  be  difficult  to  get  an  accurate  estimate  of  the  heat  generation  in  many  of  these 
substances,  even  if  a  solution  to  the  conduction  equations  were  to  be  obtained. 

3.1.2.  Convection  in  the  Resist-Substrate  Gan 

Convection  heat  transfer  caused  by  gas  in  the  exposure  station  moving  between  the  mask  and  the  resist 
could  be  a  factor  in  the  heat  transfer  away  from  the  resist.  A  simple  solution  of  the  fluid  flow  equations  will 
estimate  the  gas  velocities  in  the  gap  and  help  to  determine  the  method  and  amount  of  heat  transfer  from  the 
resist  surface. 

The  momentum  equation  is  given  by 


-Vp  +  HV^V 


(3.2) 


Neglecting  buoyancy  forces  and  simplifying  to  one  dimension  gives  the  governing  equation  for  this 


situation, 

du  d^u 

dt  dy^ 


(3.3) 


The  gas  will  be  forced  through  the  resist-substrate  gap  by  the  motion  of  the  wafer.  A  good  solution  to 
Eq.  (3.3)  would  be  obtained  by  using  a  square  wave  as  the  forcing  function;  that  is,  the  boundary  conditions 
on  the  governing  equation  would  model  the  up  and  down  motion  of  the  wafer.  However,  the  solution 
becomes  very  complicated.  An  upper  bound  on  the  velocity  can  be  found  by  assuming  the  initial  gas 
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velocity  to  be  that  of  the  wafer,  then  applying  a  negative  velocity  (in  the  opposite  direction)  for  the 
boundary  condition. 


The  nonhomogcneous  boundary  condition  at  the  walls  can  be  separated  by  splitting  the  problem  into 
steady-state  and  time-dependent  parts.  The  steady-state  solution  is  obviously  the  same  as  the  boundary 
condition  for  all  values  of  y;  viscosity  will  force  the  velocity  to  reach  a  constant  value.  Subtracting  this 
from  the  boundary  and  initial  conditions  and  solving  the  remaining  time-dependent  problem  gives  the  total 


solution: 


u{yj)=V - - ^sm;„,y  e'  '" 


where  d  is  the  gap  between  the  resist  and  the  mask,  V  is  the  velocity  of  the  mask  and  wafer,  and  the  ’s 

are  the  eigenvalues  f  ^  . 

d 

The  viscosity  may  be  calculated  from  an  empirical  relation  found  in  Keesom  (1942): 


=  5.023  (3.5) 

where  the  unit  pP  is  a  micropoise,  and  IpP  =  10*’  .  Substituting  a  temperature  of  300  K  gives  a  value 

of  201  pP .  Dividing  by  a  density  of  0.178  (Feiertag,  1994)  and  changing  units  gives  a  kinematic 
viscosity  of  1.1x10'^"'%,  about  ten  times  the  value  for  air  at  these  conditions,  due  to  an  order  of 
magnitude  difference  in  density. 

To  find  the  velocity  at  the  midpoint  between  the  mask  and  the  resist,  Eq.  (3.4)  is  solved  with  the 
values  d/  =  50pm,  y  =  25pm,  V  =  0.10”X  (Feiertag,  1994),  and  v  =  1.1  x  10"^  .  After  less  than  a 

hundredth  of  a  second,  the  velocity  is  uniform  across  the  gap.  This  means  that  there  will  be  extremely  little 
circulation  of  gas  through  the  gap,  and  convection  heat  transfer  will  be  negligible. 

Although  effects  of  the  pressure  outside  the  wafer,  including  buoyancy  effects,  might  tend  to  increase 
the  gas  velocity  in  the  gap,  it  will  be  very  small,  at  best.  In  addition,  the  oscillation  of  the  wafer  will  work 
to  keep  the  helium  from  circulating,  and  any  attempt  to  force  circulation  in  the  gap  will  jeopardize  the 
integrity  of  the  system,  as  masks  are  generally  quite  thin  and  would  be  warped  by  a  pressure  gradient  in  the 
gap. 
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The  resist  and  substrate  surfaces  are  generally  highly  polished.  The  metallic  substrate  has  a  very  high 
reflectivity;  this,  coupled  with  the  low  absorptivity  of  the  resist  (which  is  clear  PMMA),  makes  radiation  to 
and  from  the  wafer  very  small.  Because  of  the  absence  of  radiation  and  gas  movement  within  the  gap, 
conduction  will  be  the  primary  means  of  heat  transfer  from  the  resist  surface. 

The  thermal  conductivity  of  helium  in  the  irradiation  chamber  will  be  about  0.152'%.^.  In 
comparison,  the  conductivity  of  PMMA  is  about  0.19%.^  ^  ^^is  will  make  conduction  through  the  helium 
very  important.  The  direction  in  which  the  heat  is  transferred  will  depend  on  the  relative  amounts  of 
radiation  absorbed  by  the  mask  and  resist;  because  the  mask  generally  absorbs  much  of  the  radiation, 
generating  a  lot  of  heat,  the  transfer  will  most  likely  be  into  the  resist  surface. 


3.2.  Two-Dimensional  Model 

Although  wafers  are  generally  round,  the  X-ray  beam  they  are  exposed  to  is  approximately  rectangular 
in  shape.  This,  and  the  fact  that  the  beam  moves  linearly,  suggest  the  use  of  a  Cartesian  coordinate  system. 
A  solution  using  the  cylindrical  coordinate  system  would  have  to  be  three  dimensional  to  model  the 
rectangular  shape  of  the  generation  function;  this  would  result  in  a  very  complex  model.  Cartesian 
coordinates  result  in  a  relatively  simple  generation  function,  and  the  physical  system  can  be  modeled  well 
with  the  two-dimensional,  time-dependent  version  of  the  equation  for  heat  conduction  in  solids, 


d^T  ^  d^T  ^  g(x,y,t)  _  1  dT 


(3.6) 


dx  dy^  k  a  dt  '  pc 

The  thermal  conductivity  is  assumed  to  be  constant;  for  the  low  temperature  rises  expected  in  the  resist  and 
substrate,  this  is  a  valid  assumption.  The  model  is  illustrated  in  Fig.  3.1. 

Because  the  beam  is  wide  with  respect  to  its  height,  it  is  the  temperature  variations  along  its  path  (the  x 
direction  in  Fig.  3.1)  which  will  be  most  important.  A  two-dimensional  model  is  therefore  used  for  the 
wafer.  Changing  the  model  to  three  dimensions  would  add  greatly  to  its  complexity.  Because  of  the  low 
conductivity  of  the  resist  material,  most  of  the  heat  generated  by  the  X-ray  beam  in  the  resist  is  transferred 
out  through  the  substrate,  and  not  lengthwise  through  the  resist  itself.  For  this  reason,  the  boundary 
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conditions  on  the  edges  of  the  resist  matter  little,  and  using  a  two-dimensional  model  which  is  symmetric 
along  a  cross  section  of  the  beam  will  produce  the  desired  temperature  profile. 


A  Robin  boundary  condition  (boundary  condition  of  the  third  kind)  on  the  resist  surface  gives 
maximum  flexibility  in  modeling  the  heat  transfer  at  that  surface,  though  this  is  typically  used  to  model 
convection,  which  is  not  the  primary  mode  of  heal  transfer  for  the  system  studied: 

k^I^^l^+hT{x,Lj)  =  liT,  (3.7) 

dy 

Because  the  conductivity  of  the  helium  and  the  PMMA  are  comparable,  a  suitable  heat  transfer  coefficient  h 
may  be  determined  by  examining  the  conduction  resistance  across  the  gap.  Compare  the  equations  for 
conduction  and  convection  heat  transfer, 

g  =  A—  (3.8a) 

^  "'Ax 

and 


q  =  -hAAT, 


(3.8b) 
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respectively.  These  are  approximations  of  Fourier’s  law  and  Newton’s  law  of  cooling,  respectively,  where 

the  derivatives  have  been  replaced  with  finite  differences.  The  quantity  417  represents  the  temperature  drop 

across  the  gap  in  the  former  case,  and  the  difference  between  surface  and  environment  temperatures  in  the 

Ic 

second.  By  giving  h  the  value  ,  where  d  is  the  width  of  the  gap,  and  choosing  an  appropriate  AT ,  the 
conduction  situation  can  be  modeled  with  the  Robin  boundary  condition.  The  disadvantage  of  this 
approximation  is  that  the  temperatures  (and  therefore  the  heat  transfer  coefficient)  are  in  fact  functions  of 
time  and  the  heat  generation,  and  an  average  value  must  be  used. 

The  high  conductivity  and  better  cooling  of  the  substrate  (possibly  by  a  water  jacket)  justify  the  simple 
boundary  condition  for  that  surface  of 


r(x,0,r)  =  7;  (3.9) 

Making  the  assumption  of  a  good  thermal  connection  between  the  substrate  and  the  carrier,  or  support 
structure,  we  get  the  boundary  conditions 


r(0,y,r)  =  0  (3.10a) 

T{a,yj)  =  0  (3.10b) 

at  the  edges  of  the  resist.  Actually,  because  the  resist  is  so  thin,  and  because  of  the  assumption  of  no 
conduction  resistance  between  resist  and  substrate,  very  little  heat  will  be  transferred  lengthwise  through  the 
resist,  and  the  boundary  conditions  become  less  important.  Also,  the  generation  function  used  cannot  model 
the  passing  of  the  source  off  the  end  of  the  wafer,  making  these  boundary  conditions  even  less  significant. 


3.2.1.  Generation  Model 

The  heat  generation  from  the  incident  X-ray  radiation  is  approximately  an  exponential  function  of 
depth;  as  the  radiation  passes  through  the  resist,  it  loses  energy,  and  generation  decreases  (Ameel  et  al. 
1994).  This  is  illustrated  by  the  shading  in  Fig,  3.1.  Adapting  the  equations  proposed  by  Ameel  et  al. 
(1994)  to  the  present  geometry  results  in  a  generation  function  for  the  resist,  given  by 


8{x,y,t)  = 


0  Q<  x<tv 

tv<  x<  tv+  w 
0  tv+  w<  x<a 


(3.11) 
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where  v  is  velocity,  and  w  is  the  beam  width  in  the  x  direction,  as  shown  in  Fig.  3.1.  The  irradiance  Wand 
the  absorption  coefficient  fi  can  be  calculated  from  dosage  profiles  for  a  specific  application.  These 
profiles  will  depend  on  the  synchrotron  used  and  the  types  of  filters  used  to  modify  the  X-ray  radiation 
before  it  reaches  the  resist.  The  filters  can  be  tailored  to  obtain  an  optimum  generation  profile  (Ameel  et  al. 
1994). 

A  typical  dosage  profile  from  CAMD  in  Baton  Rouge  is  shown  in  Table  3.1.  Note  the  very  high  heat 
generation  in  the  silicon;  it  absorbs  far  more  of  the  X-ray  radiation.  As  noted  before,  a  variety  of  materials 
may  be  used  for  the  substrate,  and  many  have  radically  different  properties. 


Table  3.1.  Dosage  profiles  for  a  typical  wafer  at  CAMD 


Material: 

Depth  (mm) 

Dosage  (W/cm^3) 

PMMA 

0 

577 

2.0x10'^ 

489 

4.0x10'^ 

422 

6.0x10-^ 

368 

8.0x10'^ 

325 

0.10 

289 

Silicon 

0.10 

4390 

0.11 

1140 

0.12 

460 

0.13 

231 

0.15 

81 

0.20 

14 

Source:  Ameel  et  ai.  (1994). 


3.2.2.  Multiple-Pass  Solutions 

It  is  desirable  to  know  the  temperature  field  in  the  wafer  after  a  number  of  passes  of  the  X-ray  beam. 
If  a  way  could  be  found  to  modify  the  source  term  to  produce  this  effect,  it  would  probably  require  a 
different  solution  for  each  successive  pass,  and  the  complexity  would  quickly  become  prohibitive.  This 
problem  may  be  circumvented  by  including  multiple  passes  of  the  same  simple  generation  source  in  the 
mathematical  model.  Fortunately,  the  principle  of  superposition  allows  adding  of  the  temperature  fields  for 
each  pass  (Zill  1989).  The  temperature  rise  due  to  each  pass  of  the  beam  for  times  after  the  beam  has 
passed  of  the  wafer  is  represented  by  the  homogeneous  version  of  the  conduction  equation  (3.6).  That  is. 
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llie  generation  term  is  removed.  The  homogeneous  equation  is  the  basis  for  a  new  problem,  the  dissipation 
problem.  A  solution  is  found  for  one  pass  of  the  source  term,  and  used  to  determine  the  temperature 
distribution  within  the  resist  at  the  end  of  the  pass.  The  resulting  temperature  distribution  is  used  as  the 
initial  condition  for  the  dissipation  problem,  with  all  other  boundary  conditions  being  the  same.  The 
nonhomogeneous  generation  problem  can  also  be  handled  with  superposition  (Myers  1987),  and  added  to 
the  sum  of  the  temperature  fields  for  the  previous  passes.  By  adding  a  number  of  the  dissipation  solutions 
(with  times  corresponding  to  the  time  elapsed  since  their  respective  passes)  to  the  generation  solution,  the 
temperature  field  for  multiple  passes  of  the  source  is  obtained.  Because  the  model  is  symmetric,  simulating 
passes  in  opposite  directions  only  requires  reversing  the  field  for  successive  passes. 


3.3.  Two-Laver  Model 

An  analytical  problem  is  constructed  which  can  be  solved  for  the  temperature  profiles  in  two  layers 
simultaneously,  with  an  interface  condition  joining  them  at  their  common  boundary.  Although  this  two 
layer  model  is  limited  to  one  dimension,  it  should  still  be  useful  for  analysis  of  local  temperature  rises  (the 
temperature  rise  immediately  under  the  X-ray  source)  for  short  time  scales. 

For  a  two-layer  system,  temperatures  are  modeled  by  a  system  of  two  conduction  equations,  identical 
in  form,  and  given  by 


d^T, 

S\  (>')  _  1 

(97; 

k. 

dy^ 

k,  a, 

’  dt  ' 

a,  = 

Pl^l 

S2  0’)_  1 

ar. 

^2  «2 

dt  ’ 

«2=  ' 

for  0  <y<Li 

for  L,  <y<L^ 


(3.12) 


The  physical  model  is  illustrated  in  Fig.  3.2. 

To  get  the  most  flexibility  in  this  model,  convective  boundary  conditions  are  used  on  the  upper  and 
lower  surfaces: 

5r,(o,f) 


-it. 


dy 


+ /I,  r,(o,i)  =  /,,?; 
^.1  ^2(^.0  ~  ^*3  ’^1 


(3.13a) 


(3.13b) 
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where  /i,  and  /13  are  the  heat  transfer  coefficients  on  the  resist  and  substrate  surfaces,  respectively.  These 
conditions  matter  little  because  of  the  short  time  spans  this  solution  is  useful  for,  but  their  addition  adds 
little  to  the  complexity  of  the  problem.  is  the  ambient  temperature  near  the  surface  of  the  resist,  and  7) 
is  the  ambient  temperature  near  the  surface  of  the  substrate. 


0 

Li 


L2 


X-ray 

radiation 


resist 


substrate 


ik 


h2 


h3 


Fig.  3.2.  Physical  model  for  the  two  layer  problem 


At  the  interface,  the  flux  out  of  one  layer  and  into  the  next  must  be  constant,  giving  the  relation 


,  MM 


^1 


dy  dy 

The  conduction  resistance  at  the  interface  is  modeled  by  a  second  condition  on  the  interface, 

MM 


(3.14) 


-k. 


dy 


(3.15) 


This  equation  is  identical  in  form  to  a  convective  heat  transfer  boundary  condition,  but  the  constant  /12 
represents  the  inverse  of  the  conduction  resistance  across  the  interface. 

An  initial  condition  of  zero  temperature  will  be  used,  as  the  use  of  this  two  layer  model  is  limited  to 
studies  of  short  term,  localized  temperature  immediately  after  irradiation. 

Generation  of  heat  in  the  two  layers  will  be  modeled  using  the  same  exponential  representation  used 


for  the  one- layer  model,  though  there  will  be  separate  functions  for  each  layer: 
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(3.16a) 

(3.16b) 

The  irradiance  and  absorption  coefficients  are  determined  from  dosage  profiles  for  the  resist  and  substrate 
materials.  The  thickness  terms  appearing  in  the  exponents  are  needed  to  normalize  the  depth  term  and  get 
the  correct  value  of  the  exponent  for  each  generation  function. 

3.4.  Problems  with  the  Two-Laver. 

Two-Dimensional  Model 

The  conduction  equations  for  the  two-layer  problem  in  one  dimension  were  given  by  Eqs.  (3.12). 
When  these  equations  are  expanded  to  two  dimensions,  and  the  generation  term  is  ignored,  they  become 


d-T,  1  dT, 

^1 

dx‘  dy^  a,  dt 

a,  = 

Pi<^i 

d%  1  dT. 

dx^  dy^  a 2  dt  ' 

^2 

«2  = 

P2^2 

Although  the  two  layer  problem  is  well  documented  for  the  one-dimensional,  transient  case,  with 
several  explicit  solutions  in  the  literature,  the  two-dimensional  version  of  this  problem  is  not.  Ozisik  (1980) 
has  published  at  least  two  general  solutions  to  this  problem,  which  this  author  believes  to  be  incorrect. 

Ozisik  (1980)  gives  a  general  solution,  but  all  example  problems  are  for  the  one-dimensional  case.  He 
states  that  for  two  and  three-dimensional  models  the  solution  is  the  product  of  the  one-dimensional 
eigenfunctions.  These  eigenfunctions  are  obtained  by  separating  the  governing  equation  and  boundary 
conditions,  and  this  cannot  be  done  for  a  multi-dimensional  problem. 

For  the  one-dimensional  problem,  the  time  function  was  assumed  to  be  the  same  for  both  layers.  This 
allowed  separation  of  the  interfacial  conditions  and  evaluation  of  the  eigenfunctions  and  eigenvalues  for  the 
y  direction.  If  the  same  assumption  is  made  for  the  two-dimensional  problem,  the  solution  will  be  of  the 
form 


T{x,y,t)=  Xi{x)Yi{y)e{t) 


(3.18) 
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Substituting  this  into  the  governing  equations  (3.17)  gives 


a,.X,  Y^e+aiX^Yi  (3.19) 

Separation  of  variables  is  then  accomplished  by  dividing  by  the  quantity  X^Y^O,  to  give 

//  // 

Xi  Yi  6’ 

Because  j:,  y,  and  t  can  be  varied  independently,  each  of  the  terms  in  the  above  equation  must  be  equal  to  a 
constant.  Assigning  constants  to  the  two  terms  on  the  left  gives 


a, 

X: 


(3.21a) 


a,.^  =  A^ 
Y 


(3.21b) 


The  a:  direction  eigenfunctions  are  given  by  the  solution  to  Eq.  (3.21a).  It  does  not  matter  what 
boundary  conditions  are  applied;  the  general  solution  is 

^.(•*)  =  <^..m  sin  -^+  cos 4^  (3.22) 


where  X^{x)  is  the  function  from  Eq.  (3.18).  For  constant  temperature  boundary  conditions,  the 


eigenvalues  are  found  to  be 


n  fiK  f — 

Pn  = — V«/  (3.23) 

a 

Obviously,  these  eigenvalues  are  different  for  each  layer.  Regardless  of  the  boundary  conditions  used,  the 
eigenvalues  will  include  the  term  making  them  dependent  on  the  layer.  The  eigenvalues  appear  in 
the  time  term,  causing  it  to  be  layer-dependent  as  well,  thus  violating  the  assumption  made  at  the  beginning 
of  the  derivations.  The  practical  effect  of  this  is  that  the  temperatures  7;  and  do  not  agree  with  the 
boundary  conditions  at  the  interface.  For  instance,  the  temperature  might  have  an  unusually  large 
discontinuity  for  very  low  interfacial  resistances,  or  might  jump  in  the  wrong  direction,  giving  a  temperature 
profile  that  is  physically  impossible.  This  is  how  the  problem  was  discovered,  after  an  (incorrect)  solution 
had  been  obtained  and  the  evaluation  program  debugged. 


20 


Dropping  the  assumption  of  identical  time  terms  causes  problems.  The  interfacial  conditions  are  no 
longer  separable,  so  individual  problems  for  the  A:(jr)  and  y(>’)  eigenfunctions  cannot  be  created.  Not 
being  able  to  find  these  eigenfunctions  precludes  their  being  multiplied  to  obtain  a  solution,  as  Ozisik 
instructs  the  reader  to  do  in  Heat  Conduction. 

In  a  telephone  conversation,  Ozisik  (1995)  admitted  that  the  problem  was  “more  complicated”  than  his 
text  implied.  He  said  he  was  aware  of  only  a  few  researchers  who  had  succeeded  in  solving  multiple-region 
problems  for  more  than  one  dimension,  and  that  each  of  these  solutions  required  major  simplifications  to 
the  problem.  One  example  of  this  type  of  solution  is  a  paper  written  by  Sareen  and  Gidaspow  (1974)  on  the 
subject  of  mass  diffusion.  The  model  used  in  the  paper  is  a  two-layer,  two-dimensional  slab  much  like  the 
one  originally  envisioned  for  the  LIGA  problem.  This  looked  like  a  promising  breakthrough  until  it  was 
realized  that  the  author  had  assumed  the  diffusion  coefficient  to  be  the  same  in  the  lengthwise  direction  for 
both  layers.  For  the  problem  outlined  above,  this  is  the  equivalent  of  letting  the  thermal  diffusivity,  a,  have 
the  same  value  for  both  resist  and  substrate  in  the  x  direction  (a  very  bad  assumption  for  the  LIGA 
situation).  Assuming  equal  diffusion  coefficients  in  the  two  layers  causes  the  jr-direction  eigenvalues  to  be 
the  same  for  both  layers,  which  in  turn  makes  the  time  terms  the  same,  and  the  initial  assumption  is  not 
invalidated.  This  peculiarity  allowed  the  problem  to  go  undetected  for  a  long  time.  Among  the  arbitrary 
values  used  to  test  the  (incorrect)  solution  obtained  from  Ozisik’s  general  solution  were  identical  values  of 
thermal  diffusivity,  even  while  it  was  being  tested  with  different  values  of  thermal  conductivity,  which  does 
not  cause  a  problem.  It  was  only  when  actual  material  properties  were  used  and  a  number  of  cases  studied 
that  the  problem  was  recognized. 


CHAPTER  4 


MATHEMATICAL  SOLUTION 

This  chapter  will  describe  the  solution  of  the  mathematical  models  of  the  resist-substrate  system.  The 
validity  of  the  solutions  will  be  checked  by  comparing  them  with  simpler  solutions,  and  by  an  analysis  of 
the  data  they  produce.  For  instance,  a  steady-state  solution  may  be  used  to  confirm  the  results  of  a  time- 
dependent  solution  by  solving  the  latter  for  large  values  of  time. 

4.1.  One-Laver  Solution 

Because  of  the  form  of  the  generation  term,  Eq.  (3.10),  only  one  pass  of  the  X-ray  source  can  be 
modeled  with  any  one  solution.  It  will  be  desirable  to  know  the  temperature  distribution  in  the  resist  after  a 
number  of  passes,  and  a  single  generation  term  which  sweeps  back  and  forth  over  the  resist  would  be 
extremely  complicated.  The  principle  of  superposition  can  be  used  to  add  the  effects  of  individual  passes  of 
the  generation  term  to  obtain  a  final  temperature  distribution.  To  find  the  temperature  during  the  nth  pass, 
solutions  are  found  not  only  for  that  pass  but  for  all  the  ones  that  preceded  it.  Only  the  present  pass  has  a 
generation  term;  all  others  arc  simply  time-decay  solutions,  with  an  initial  condition  consisting  of  the 
temperature  distribution  at  the  end  of  one  pass  of  the  source  term. 

The  final  solution  for  multiple  passes  of  the  X-ray  source  will  utilize  two  similar  versions  of  the 
problem.  The  first  will  include  the  generation  term  and  will  be  used  to  calculate  temperatures  for  a  single 
pass  (the  present  pass)  of  the  X-ray  source.  The  initial  condition  for  this  solution  will  be  zero  temperature 
in  the  wafer.  The  second  version  of  the  solution  will  be  used  to  calculate  the  decay,  or  cooling,  of 
temperature  from  previous  passes  of  the  source  over  the  wafer.  The  initial  condition  for  this  version  will  be 
the  temperature  profile  after  one  complete  pass  of  the  source,  and  the  problem  will  have  no  generation  term. 
To  obtain  the  temperature  distribution  after  a  number  of  passes,  the  distribution  from  the  generation 
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problem  will  be  superimposed  on  a  series  of  decay  solutions  solved  for  an  appropriate  time,  corresponding 
to  the  time  elapsed  since  each  of  the  previous  passes.  Alternating  the  direction  of  the  decay  solution  will 
model  the  passing  of  the  X-ray  source  in  opposite  directions;  because  the  resist  is  symmetrical,  the  same 
solution  can  be  used  for  both  directions,  and  the  grid  of  solution  data  can  simply  be  reversed. 

The  solution  of  the  conduction  equation  in  one  layer  is  straightforward,  and  well  documented  in  many 
conduction  texts  (Ozisik,  1980  and  Myers,  1987). 

The  complete  statement  of  the  one-layer  problem  with  general  generation  and  initial  condition  terms  is 

d'^T  g(x,y,t)  I  k 

(4.1) 

with  the  boundary  conditions 

r(x,0,r)  =  0  (4  2a) 

dT{x,L,t) 

—  +  HT(x,L,t)  =  0  (4.2b) 

T{0,y,t)  =  0  (4  2c) 

T[a,y,t)  =  0  (4.2d) 

where  H  is  the  heat  transfer  coefficient  h  divided  by  the  conductivity  K  and  the  ambient  temperature  is 
assumed  to  be  zero.  Thus,  T  may  represent  a  temperature  rise  in  the  resist  above  . 

The  initial  condition  for  a  single  pass  of  the  generation  term  is 

Ti{x,ySi)  =  0  (4  3) 

and  the  generation  term  is 


g{x,y,t) 


0  0<  jc  <  rv 

tv^x<tv+w 
0  tv+  w<  a 


For  the  time-decay  solutions,  the  generation  term  is  g{x,  y,  f)  =  0 ,  and  the  initial  condition  is 


(4.4) 


V  V  J 


(4.5) 
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where  the  quantity  -  is  the  time  needed  for  the  leading  edge  of  the  X-ray  beam  to  reach  the  end  of  the 
wafer;  the  beam  is  assumed  to  travel  the  length  of  the  wafer,  minus  its  own  width.  If  t  is  allowed  to  become 
large  enough  that  part  of  the  beam  could  move  beyond  the  outer  edge  of  the  wafer,  the  solution  is  incorrect; 
this  condition  does  not  correspond  to  the  physical  situation. 

A  separation  of  variables  solution  of  Eq.  (4.1)  is  begun  by  assuming  that  the  solution  is  of  the  form 


7’(ac,y,/)  =  K(y)X(A-)0(r) 

(4.6) 

Substitution  into  the  governing  equation  (4.1)  yields 

x're+xY"e=—XYe' 

a 

(4.7) 

Dividing  this  equation  by  XY6  gives 

X'^  Y"  _  \  6' 

X  ^  Y  “aT 

(4.8) 

For  Eq.  (4.8)  to  hold,  each  of  the  terms  must  be  equal  to  a  constant. 

This  allows  us  to  separate  the 

equation  into  the  three  problems: 

X  ^ 

(4.9a) 

Y" 

—  =  -A^ 

Y 

(4.9b) 

a  e 

(4.9c) 

The  signs  of  the  constants  in  Eqs.  (4.9)  are  interchangeable,  as  are  the  forms  of  the  general  solutions 
derived  from  these  equations.  The  choice  is  not  arbitrary;  only  one  will  lead  to  the  correct  solution.  Here, 
the  choice  is  made  with  a  previous  knowledge  of  which  solution  will  lead  to  an  answer. 


The  general  solution  to  the  X(j:)  part  of  the  problem  is 

X{x)  =  Acos  fix +  B sin  Px  (4.10) 

By  substituting  Eq.  (4.6)  into  the  boundary  conditions  given  by  Eqs.  (4.2c  and  d)  and  dividing  by  the 
unchanged  terms,  the  boundary  conditions  on  Eq.  (4.10)  arc  obtained  as 
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A'(0)  =  0 

X(a)  =  0 


(4.11) 


Applying  the  first  of  these  boundary  conditions  to  Eq.  (4. 10)  shows  that  A  is  equal  to  zero.  The  second  one 
gives  the  relationship 


B  sin  Pa  =  0 

For  this  equation  to  hold,  the  eigenvalues  P  must  be: 


p„=—,  n  =  0,1,2,. 
a 


and  the  eigenfunction  for  the  x  direction  becomes 


(4.12) 


(4.13) 


X{x)=Bsin^^^  (4.14) 

a 

The  general  solution  of  the  separated  equation  in  y  (4.9b)  is 

K(v)  =  Csin  A>’4*  DcosA>'  (4T5) 

The  constants  C  and  D  are  determined  by  applying  the  boundary  conditions  in  the  y  direction. 
Equation  (4.6)  is  substituted  into  the  boundary  conditions  (4.2a  and  b),  resulting  in  the  boundary  conditions 
for  the  function  T(y) : 

y(0)  =  0  (4.16) 

^^^^+HY{L)  =  0  (4.17) 

dy 

Applying  the  first  of  these  conditions  to  the  general  solution  (4.15)  results  in  the  condition  D  =  0. 
Application  of  the  second  condition  gives 

A  cos  AL  +//  sin  AL  =  0  (4. 1 8) 

For  this  equation  to  hold,  the  eigenvalues  A^  must  be  the  roots  of  the  equation 

-A^cot  A^L=//  (4.19) 


The  y  direction  eigenfunction  is  then 


f(y)  =  CsinA„y 


(4.20) 


By  substituting  the  separated  solutions  (4.20  and  4.14)  into  Eq.  (4.6),  the  final  solution  takes  the  form 
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n.x,y,t)  =  J^'^A(t)sinli„xs\nX,„y  (4.21) 

n=l  m=l 

with  and  defined  by  Eqs.  (4.13  and  4.19),  respectively.  The  constants  A  and  C  have  been  integrated 
into  the  function  44(/) . 

The  above  solution  is  substituted  into  Eq.  (4.1),  and  orthogonality  is  applied.  In  the  jc  direction,  the 
orthogonality  relation  is 


f"  .  nnx  ,  mnx  , 

sm - sin - dx- 

Jo  a 


m^n 

m 


The  orthogonality  relation  for  the  y  direction  is 


(4.22) 


.  ;i;ry  .  m;ry  , 
sin - sin - —dy  = 


l[x\+h^)+h^ 


m^n 


m=n:itO 


(4.23) 


The  simplifications  used  to  obtain  the  quantity  on  the  right  side  of  this  equation  utilized  Eqs.  (4.17  and 
4.15),  and  are  relatively  simple,  if  not  immediately  obvious.  Ozisik  (1980)  presents  an  example  of  this 
simplification  for  a  general  solution  with  the  boundary  conditions  used  here. 

The  orthogonality  constants  given  by  Eqs.  (4.22  and  4.23)  may  be  combined  into  a  single  constant  for 
this  solution,  N„„: 

’  n,m 


N„.„.  = 


4{Xi+H^) 


(4.24) 


The  general  solution  (4.21)  is  now  substituted  into  the  governing  equation  (4.1).  Orthogonality  is  applied 
by  multiplying  by  the  appropriate  eigenfunctions  and  integrating,  to  get 


-A{t)PlNn_„  -  A{t)X‘„N„^  +  £  £  iiili:llsin  A„y  sin  li„x  dydx  = 


(4.25) 


Rearranging: 


CC  [i-gix,y,t)  a  .  . 

- - sin  A„ysin  P„x  dydx 


(4.26) 
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The  integration  on  the  right  hand  side  eliminates  the  dependence  on  x  and  y,  leaving  a  nonhomogeneous 
differential  equation  which  may  be  easily  solved  for  the  function  A{t) .  By  examination,  the  integrating 

factor  is  observed  to  be  .  Multiplying  by  this  factor  and  integrating  both  sides  with  respect  to  t 

gives  the  relation 


a  f  [•“  r'-  «(-r.  J’.o 

Nn..  J  Jo  Jo  k 


sin  A„,  y  sin  dydxdt 


(4.27) 


and  the  solution  to  the  time  dependent  part  of  the  problem  is  given  by  Eq.  (4.21)  with  A{t)  given  by  the 
relation 


^10  ^  e  je  ^  J  - - - sin  X,„y sin  p„x  dydxdt 


(4.28) 


where  /l(/)  is  seen  to  be  dependent  on  the  general  function  g(x,  y,r) ,  and  the  initial  condition  will  be  used 


to  complete  the  indefinite  integral  with  respect  to  time. 

4.1.1.  First  Pass  Solution 

To  obtain  the  solution  for  a  single  pass  of  the  X-ray  source  over  the  resist,  Eq.  (4.28)  is  solved  with 
the  initial  condition  and  generation  functions  given  by  Eqs.  (4.3  and  4.4),  respectively.  The  integral 
involving  the  generation  term  is  solved  first. 

The  generation  term,  Eq.  (4.4),  specifies  that  the  generation  is  zero  everywhere  but  on  the  interval 
tv<x<tv+w.  Therefore,  the  integral  with  respect  to  x  in  Eq.  (4.28)  is  zero  for  all  x  not  on  that  interval, 
and  the  integral  reduces  to 


[“  g{x,y,t)  .  „  ,  r"'+»'5(x,y,/)  .  „ 

Jo-^sin^„x^  =  J^  ^-^sinP„x 


g{x,y,t) 


k  ’  “  Jfv  ic 

The  X  dependence  of  the  generation  function  has  been  removed,  and  the  int 


(4.29) 


r+>vg(x,y,t)  .  „  ,  gUy.t)  -Ir  „  / 

— - sinP„xdx=-^  ^  [cos^„(tvf 


The  integral  with  respect  to  y  is  straightforward: 


.5) 
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rLWfte^''-' 

Jo - 1 - sinA„)'d^ 

Wn  ,  •  / 

k  4,2+ 32  A„cosA„>-) 

M  +  >^m  Q 

W^  \  /  ^ 

k  fl^  •¥  }}  ^ 


-(/zsinA^y-  A„cosA„>-) 


Substituting  these  results  back  into  Eq.  (4.28),  we  get 


4)  / 1  (■«(/»;+ Ai) ,  v^__^2__ 

(/isinA„L-  A„cosA„L+e  ^^A,„j[cos)3„(fv  + w)-  cos/?„fv]c//| 


and  performing  the  integration  with  respect  to  time  gives 


..'1.  .. .  <«3) 


^nm  k  /J(u-  +  AM''  '"  "'  la2  ,2^2  nil 

"■"L  H^.+A^jj  +^y 

[a[Pl  +  A^„  )[cos  P„{tv  +  w)-  cos  p„tv]  +  )3„  v[sin  P„(tv+w)- sin  /)„fv]|  -  e'  ^"'^'  F{x,  y) 
where  F{x,y)  is  a  “constant”  of  integration. 

Applying  the  initial  condition  (4.3)  to  Eq.  (4.21)  leads  to  the  requirement  that  the  function  A(r)  equal 
zero  for  t  =  0;  this  is  the  only  way  the  summation  can  be  zero  for  all  points.  The  function  F{x,  y)  is  then 
easily  determined: 


cosA^L+e-'-Uj 

H/i; + A- )]■  *  A' 


and  /\(/)  becomes 


fi  f„2  .  22  A„cosA„L+e-'"-A„)-— - —5 - — 

|a(  /J^  +  A^,„  jj^cos  )3„  (/V + w)  -  cos  /J„/v  +  e"  ^ '  (cos  /3„  w  - 1 )  j  + 

/J„ vj^sin  P„ {tv  +  w)  -sin  p„tv + e'  sin j3„ w  | 
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This  completes  the  solution  for  one  pass  of  the  X-ray  source  over  the  resist.  Verification  of  this 
solution  will  be  accomplished  using  the  two-layer  solution  at  the  end  of  this  chapter. 

Figure  4.1  shows  a  temperature  distribution  produced  by  the  one  pass  solution,  using  the  constants  in 
Table  4.1. 


5.4 

Fig  4.1.  Temperature  rise  for  generation  part  of  one-layer  solution 


Table  4.1.  Constants  used  for  generation  solution 


L 

0.03  cm 

a 

6.0  cm 

h 

100 

nr  K 

a 

1.182x10'’ 

k 

0.198  * 

W 

250360 

nr 

V 

0.10  f 

2304.7  iL 

t 

3.0  s 

w 

1.0  cm 

These  constants  were  chosen  for  illustrative  purposes,  but  most  correspond  to  the  conditions  one  might  find 
in  an  exposure  chamber.  The  material  properties  a  and  k  and  the  thickness  L  are  typical  for  a  PMMA  resist. 
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Likewise,  the  generation  function  constants  W  and  /t  correspond  to  a  filtered  exposure  at  CAMD  with  a 
synchrotron  energy  of  1.5  GeV  and  moderate  current.  The  length  of  the  resist  and  the  width  of  the  X-ray 
beam  are  both  unrealistically  long;  these  values  make  the  trend  of  the  temperatures  in  Fig.  4.1  more 
obvious. 

In  Fig.  4. 1 ,  the  X-ray  beam  is  moving  from  left  to  right.  The  temperature  is  zero  on  the  right  half  of 
the  graph,  where  the  resist  has  not  yet  been  exposed,  and  it  rises  as  the  beam  passes  over  the  center.  The 
temperature  is  highest  at  the  trailing  edge  of  the  beam,  which  has  received  the  maximum  exposure,  then 
drops  on  the  left  side  of  the  graph  due  to  conduction  and  convection  out  of  the  resist.  The  steeper  slope  on 
the  left  end  of  the  graph  is  caused  by  the  beam  having  started  there;  a  full  dose  of  radiation  is  not  received 
for  those  points.  The  temperature  at  the  back  of  the  graph  drops  to  zero,  as  required  by  the  boundary 
condition  given  by  Eq.  (4.2a). 

4.1.2.  Time-Decay  Solution 

The  first  step  in  the  solution  of  the  time-decay  problem  is  to  determine  the  temperature  distribution  at 
the  end  of  one  pass  of  the  source  term.  The  appropriate  length  of  time  for  the  end  of  the  first  pass  is 
^  ^  prime  is  added  to  the  time  variable  in  subsequent  formulas  to  distinguish  it  from  the  time 

used  in  the  most  recent  pass,  which  is  used  in  the  generation  solution.  However,  when  the  complete, 
superimposed  solution  is  evaluated  for  a  number  of  passes,  the  decay  solution  will  be  solved  with  different 
times  corresponding  to  each  of  the  passes.  Letting 


will  simplify  later  derivations.  This  will  make  the  time  elapsed  since  the  end  of  the  pass  whose 
generation  the  solution  represents.  Substituting  Eq.  (4.36)  into  Eq.  (4.21),  with  A(r)  given  by  Eq.  (4.35), 
gives 

TfAT,  y,  /'  =  0)  =  ]£  ]£  a[  sin  sin  A„y  (4.37) 

nsl  m=l  V  V  y 


where 
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4(r-=0)  =  ^fe.  (//sin  A„,L- cosA„,L  +  .-''"A„,) 


+  cosP„a-cosp„{a-w)  +  e 


(cosi3„w-l) 


P„v  s\nP„a-s\nP„{a-w)  +  e  »  sin^„w 


Because  there  is  no  generation  for  this  part  of  the  solution,  g(j:,  y,t)  =  0,  and  the  general  coefficient 
given  by  Eq.  (4.28)  simplifies  to 


where  the  function  B[t')  will  be  the  series  coefficient  for  the  time-decay  solution.  Equation  (4.38)  is  the 
initial  condition,  and  is  used  to  solve  for  the  function  Fj[x,y)  by  equating  it  with  the  time-decay  solution, 
with  r  =  0.  Because  both  functions  arc  double  series  of  the  same  form,  the  coefficients  must  be  equal;  this 


requirement  gives 


V  V 


Solving  for  substituting  into  Eq.  (4.39),  we  get 


and  the  solution  to  the  decay  part  of  the  problem  is 


=  sin)3„j:sinA„y 


The  above  equation  is  exactly  what  one  would  expect  for  this  solution;  intuitively,  it  makes  sense.  The 
time-decay  solution  is  simply  an  exponential  decay  of  the  temperature  distribution  when  the  heat  source  is 
taken  out  of  the  problem  (corresponding  to  its  passing  off  the  end  of  the  wafer). 

A  plot  of  the  temperature  as  a  function  of  time  is  shown  in  Fig.  4.2.  The  variables  used  in  this 


calculation  are  from  Table  4.1,  with  the  exception  of  time,  which  has  been  changed  to  60  seconds. 
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Fig,  4.2.  Temperature  for  multiple  passes  of  the  X-ray  beam  at  center  of  resist 


Each  of  the  small  oscillations  of  the  temperatures  in  Fig.  4.2  corresponds  to  one  pass  of  the  X-ray 
beam.  During  the  initial  passes,  the  temperature  increases  quickly  from  one  pass  to  the  next.  After  60 
seconds,  this  increase  has  slowed,  and  the  temperatures  approach  a  “steady-state”  oscillation  pattern,  with 
constant  temperature  range  and  maximum  temperature. 


4.2.  Steady-State  Solution 


If  multiple  solutions  are  not  going  to  be  added  together,  nonhomogeneous  boundary  conditions  can  be 


considered.  This  may  be  used  to  compare  the  effects  of  different  types  and  amounts  of  heat  transfer  at  the 
surfaces  of  the  resist,  and  the  flow  of  heat  through  and  along  the  resist.  The  steady-state  problem  in  two 
dimensions  is  stated  as 


5  V  ^  V 
^  =  0 
dx^  dy^ 

with  the  boundary  conditions 


(4.43) 


(4.44a) 


d^{x,L) 

dy 


+  H'F{x,L)  =  HT^ 


(4.44b) 
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'P(0,y)  =  0 

(4.44c) 

f'(«,y)  =  0 

(4.44d) 

Separation  of  variables  is  begun  by  assuming  a  solution  of  the  form 

'I'{x,y)  =  Y{y)X{x) 

(4.45) 

Substituting  Eq.  (4.45)  into  the  governing  equation  (4.43)  and  dividing,  we  get 

X"  Y" 

(4.46) 

Each  of  the  terms  on  the  left  side  of  this  ec^uation  must  be  a  constant,  because  x  can  be  varied  independently 
of  y,  so  we  get  the  two  equations 


(4.47a) 

Y"  2 

—  =  1 

Y 

(4.47b) 

with  the  general  solutions 

X(a')  =  Asin  rjx  +  Bcos  rjx 

(4.48a) 

K(y)  =  Csinh  ryy  +  Dcosh  Tfy 

(4.48b) 

Applying  the  boundary  conditions  in  the  x  direction  shows  r]  to  be  ,  where  n 

a 

is  the  set  of  positive 

integers.  Also,  the  coefficient  B  must  be  zero.  The  A-direction  eigenfunction  becomes 

X(a)  =  Asin77„A 

(4.49) 

multiplying  the  solutions  given  by  Eqs.  (4.49  and  4.48b),  and  letting  AC  =  C„ ,  and  AD  =  D„,  we  get 

(^>  y)  =  S s>"  sinh  Ti^y  +  D„  cosh  i]„ y)  (4.50) 

«=! 

applying  the  boundary  condition  at  y  =  0  gives 

'f' =  sin  ri„x  =  T, 

n=l 


(4.51) 
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This  is  in  the  form  of  a  Fourier  series  representation  of  T,,  which  allows  the  constant  D„  to  be  easily  found. 
From  Greenberg  (1988),  we  get 


=-[jismr]„xdx 
a  '^0 


(4.52a) 


Applying  the  convection  boundary  condition  (4.44b)  to  Eq.  (4.50)  results  in  the  solution  for  C  : 


(4.52b) 


^sinri„x{C„  r/„  cosh77„L+  D„  t]„  sinh77„L)  + 


(4.53a) 


V„x{C„  sinh  T]„L  +  D„  cosh  7)„L)  =  HT^ 


Rearranging  and  again  using  the  Fourier  representation,  we  get 


cosh t]„L+H sinh  r]„L)  +  sinh 77„L  +  W cosh  ^  sin  ri„x  dx  (4.53b) 


and  substituting  the  solution  for  given  by  Eq.  (4.52b)  results  in: 


~  1  “ (-  0" ][^^u  -  T; (i7„  sinh r]„L+H cosh t]„ L)] 

=  - - - — -  (A 

77„  cosh  77„L+// sinh  7]„L  ^  ‘  ^ 

The  solution  to  the  steady  state  problem  is  given  by  Eq.  (4.50)  with  the  constants  given  in  Eqs.  (4.52b 


and  4.54). 


4.2.1.  Confirmation  of  the  Steady-State  Solution 


To  check  the  steady-state  solution,  a  one-dimensional  problem  is  solved,  with  boundary  conditions 
identical  to  the  y-direction  boundary  conditions  on  the  two-dimensional  problem.  The  governing  equation 


is  then 


with  the  boundary  conditions 


y{o)=T, 


(4.56a) 


dY{L) 


HY{L)  =  HT, 


(4.56b) 
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Equation  (4.55)  is  integrated  twice  with  respect  to  j,  and  the  solution  takes  the  form  of  the  linear  equation 


Y{y)^Ay-^B 


(4,57) 


Applying  the  boundary  conditions  and  solving  for  the  constants  A  and  B  gives 


A  = 


l  +  L// 


(4.58a) 


5=7; 


(4.58b) 


Observe  that  the  two-dimensional  solution,  Eq.  (4.50),  must  agree  with  the  boundary  conditions  in  the 
X  direction,  because  each  of  the  sin  T]^x  terms  in  the  summation  is  zero  when  x  is  equal  to  0  or  a.  The 


one 


dimensional  solution  above  is  used  to  check  agreement  with  the  remaining  boundary  conditions  in  the  y 
direction,  and  the  upper  and  lower  surfaces  of  the  resist.  The  resist  is  made  to  approximate  an  infinite  slab 
by  making  it  far  longer  than  it  is  thick.  Points  near  the  middle  will  be  affected  very  little  by  the  heat  loss 
from  the  ends  of  the  resist,  and  the  temperature  profile  of  a  cross  section  should  be  very  close  to  the  one¬ 
dimensional  solution.  Comparison  of  the  two  solutions  will  also  check  the  computer  program  used  to 
evaluate  them. 

Table  4.2  shows  the  properties  and  dimensions  used  for  the  calculations.  The  boundary  temperatures 
are  arbitrary,  but  choosing  them  to  be  nonzero  minimizes  the  chance  of  missing  errors  in  the  solutions.  The 
long,  thin  resist  (a  and  L)  causes  the  two-dimensional  solution  to  approximate  the  one-dimensional  solution. 


Table  4.2.  Constants  for  test  of  the  steady- 
state  solution 


L  0.01  cm 

t;  5k 

a  6.0  cm 

T,  2K 

k  0.198  * 

h  1600  -IV 

nr  K 

The  one  and  two-dimensional  solutions  both  produce  a  temperature  profile  in  the  form  of  a  straight 
line.  Table  4.3  lists  the  values  of  the  two  functions  at  the  top  and  bottom  surfaces,  showing  very  good 
agreement  at  these  boundaries.  The  discrepancies  in  the  two  values  are  due  to  conduction  out  of  the  sides 
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of  the  two-dimensional  slab  (though  this  is  very  small),  and  truncation  of  the  infinite  sum.  The  latter  topic 
will  be  discussed  in  more  detail  in  Chapter  5. 


Table  4.3.  Boundary  point  temperatures  for  steady-state  solution 


Surface: 

Temperature  Rise 

2-D  solution 

1-D  solution 

top 

3.34043 

3.34078 

bottom 

1.99886 

2.00000 

The  temperatures  produced  by  the  two-dimensional  solution  vary  little  in  the  x  direction,  because  of 
the  extreme  thinness  of  this  resist.  However,  they  do  decrease  toward  the  outer  edges  of  the  resist,  and  drop 
to  zero  at  the  boundaries,  as  expected. 


4.3.  Two-Layer.  Time-Dependent  Solution 

The  one-dimensional  problem  in  two  layers,  outlined  in  Chapter  3,  is  stated  as 


<5=7;  ^ 

S\  iy)  _  1 

dT, 

for  0<y<Li 

k,  or, 

'  dt  ' 

a,  = 

Pi^i 

(4.59a) 

I 

dT^ 

for  L, 

dy-  ^ 

^2 

dt  ' 

«2  =— ^ 
P2C2 

(4.59b) 

with  the  boundary  conditions 


dTJ0,t)  ,  ^ 

(4.60a) 

(4.60b) 

^dr,(L„t)  dr2m,t} 

dy  dy 

(4.60c) 

I  7/  7  f  r  A  77  T' 

(4.60d) 

the  initial  condition 
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7;(.v,0)  =  0 


and  the  generation  terms 


(4.62a) 


Mr-U-k) 


(4.62b) 


S2{y)='Wili2e 


The  following  simplifications  have  been  made  for  conciseness; 


//,  =-5^;  H,=^-  K=^ 


To  begin,  it  is  assumed  that  a  solution  exists  of  the  form 

7:(3'.0=>:0')^(0  (4.64) 

Note  that  the  time  function  9{t)  is  the  same  for  both  layers.  Substituting  this  into  the  homogeneous  form  of 
the  governing  equations  (4.59a  and  b)  and  dividing  by  Tj(y)0(f)  gives 

// 

Y;  6' 

«/-]r  =  7  (4.65) 

The  two  sides  of  this  equation  must  be  equal  to  a  constant.  Because  the  time  function  was  assumed  to  be 
the  same  for  both  layers,  the  terms  on  the  left  must  be  equal  to  the  same  constant,  though  their  component 
functions  are  layer  dependent.  This  requirement  gives 


Y 


The  general  solution  to  this  equation  is 


My)  =  sin  S,>,  cos^^ 


Note  that  the  above  is  not  the  same  as  the  one  used  in  the  two-dimensional  solution;  in  fact,  it  has 
different  dimensions.  The  eigenvalues  ( )  and  corresponding  specific  eigenfunctions  (  ^  and  5.  ^ )  are 

found  by  substituting  the  general  form  of  the  eigenfunction  given  by  Eq.  (4.67)  into  the  homogeneous  form 
of  the  boundary  conditions  (4.60),  resulting  in  the  system  of  equations: 
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^l,m  I - 

V«1 


(4.68a) 


A  R 

^\,m  i - COS— -  Di 

a/«,  Ja, 


■ 

— sin 

./I  d 

^2,m  I -  “  "2,m  COS  — 


4  A,^4^cos^-5.„.4^sin^  U  /i. 


a,  Ja 


a, 


^fn  — n  -^m^l 


CCl  ^cc 


ccj  Ja 


(4.68b) 


(4.68c) 


cosi^-  B2.„  -^s\n^^+  h\  sin^^+  cos 

^CCi  ./cK-i 


=  0  (4.68d) 


This  set  of  equations  can  be  expressed  more  concisely  in  matrix  form.  Using  the  simplifications 


7  =  and  7]  =  -~ 

V«i  Ja 


and  grouping  the  terms  results  in 


.  7cos)L,  ^  7sin)L, 
sm  tL,  +4 - LJ.  cos  ~4 - L± 


^7  cos  yL^ 
0 


-^7 sin 
0 


0 

-sin  rjL, 

-7]  cos  r\L^ 
77COST7L2 


0 

-cos  rjL^ 
77  sin  rjL^ 


0 


OS77L2  .  ,  T^sinT]!, 

- - +sin77L2 - -^+cos7]L,  0 

“3  4/3  L  •  J  U  J 


Because  these  equations  are  homogeneous,  the  constants  A^  „ ,  ^ ,  B,  „ ,  and  Bj  „  can  only  be  found  in 

terms  of  any  one  of  the  four,  or  in  terms  of  some  arbitrary  constant.  However,  this  is  not  a  problem, 
because  in  the  final  solution,  they  will  appear  in  both  the  numerator  and  denominator,  and  the  arbitrary 
constant  will  cancel  out.  With  this  in  mind,  then,  the  constant  A]  is  set  equal  to  one,  and  Eqs.  (4.68b,  c, 
and  d)  are  used  to  find  the  remaining  constants: 


=  1 


(4.71a) 


(4.71b) 
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^2,m  - 


TysinrjL,  +511177^2  |+  7? cos t?L,  -  +COS77L2 

V  "3  y  I  W3 


(4.71c) 


Ba,«  = 


I'J  rsin>i.,  Y  ncosTjz,  ,  1 

Ay  COS)Lj  —  hSinTjZ/j 

_ V _ ^1  A  “3 _ / _ 

77sin7)Z,|  — +sinriL^  +r|cosr|L^  -  +cost)L; 

^  "3  y'  I  //j 


(4.71d) 


The  equation  for  the  eigenvalues  is  obtained  from  the  requirement  that  the  determinant  of  the 
coefficient  matrix  in  Eq.  (4.70)  must  be  zero,  as  shown  in  Eq.  (4.72). 


si„,t,+Ii22A  cos,t,-Iiii4i 

//-)  H-j 


Ky  cos  yL^ 


-Ky  sin  yL^ 


‘Sin  7]L, 


-Tjcos  tjL, 


-  cos  rfL^ 


7]  sin  77L, 


=  0  (4.72) 


rfcosrfLj  .  ,  risinriLy 

— -^+sin7]L2  ^  +COS77L, 

“3  “3 


This  completes  the  solution  of  the  eigenvalues,  given  by  the  roots  of  the  transcendental  equation 
(4.72),  and  the  eigenfunctions,  given  by  Eq.  (4.67)  with  the  constants  in  Eqs,  (4.71a,b,c,  and  d).  The 
boundary  conditions  are  now  used  with  these  eigenfunctions  to  obtain  a  complete  solution. 

Ozisik  (1980)  obtained  the  general  solution  to  a  two-layer  problem  with  generation  using  integral 


transforms: 


(4.74a) 


(4.74b) 


(4.74c) 


In  this  notation,  F.(r')  is  the  initial  temperature  distribution  in  layer;,  /;(r',r')is  the  temperature 
distribution  on  the  surface  of  layer;,  and  the  general  dimensional  variable  r  is  used.  In  the  present  solution 
this  IS  replaced  by  the  single  thickness  variable  y.  For  the  resist  surface,  /|(r',r')  =  //,?; ,  and  for  the 
substrate  surface,  f2{r',t')=  H^Tf.  The  initial  condition  (4.61)  states  that  the  initial  temperature 
distribution  is  zero  for  the  whole  volume.  This  makes  the  term  F(A,„)  zero,  and  Eqs.  (4.73  and  4.74) 


reduce  to 


where 


fn=\ 

j-\ 


(4.76a) 


(4.76b) 


The  orthogonality  constant  reduces  to 

and  substituting  the  appropriate  limits  for  the  present  problem  gives 

=  -—r'  ^ y^Ayfdy 

For  the  eigenfunctions  under  consideration,  the  integral 


is  given  by 


J  dy 

I  Va,  Ja, 


Integration  of  this  term  gives 


aJ  ^2  +  (4.81) 


Substituting  the  above  equation  into  the  definite  integral  in  Eq.  (4.78)  and  using  the  simplifications  given  in 
Eq.  (4.69)  gives  the  final  form  of  the  orthogonality  constant: 
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k,  n 


"■^{[2  (^'•'"  ■"  -  0^‘"2;7Z.  -  COS^  77Z,  1 


-('^J.-n  +^2,m))^i  +— (52,„  -  Aj^)sin2}f,|  cos^  yL^ 


The  boundary  condition  terms  in  Eq.  (4.76)  are  easily  evaluated 


(4.83a) 


-  ^3^(^2.™  sinrjL^  +  ,„  cosrjZ^)  (4.83b) 

These  terms  are  in  an  integral  with  respect  to  time  in  Eq.  (4.75).  Because  they  do  not  vary  with 
respect  to  time,  they  may  be  removed  from  this  integral,  which  is  then  evaluated  as: 


At  this  point,  enough  of  the  solution  has  been  constructed  to  allow  a  simple  test  of  the  derivations. 
Using  the  boundary  condition  terms  in  Eqs.  (4.83a  and  b),  the  orthogonality  function  in  Eq.  (4.82),  and  the 
time  decay  term  (4.84),  Eq.  (4.75)  can  be  evaluated  to  get  some  actual  temperature  rises. 

4.3.1.  Test  of  Solution  With  No  Generation 

Although  the  two-layer  solution  is  less  calculation  intensive  than  the  multiple-pass,  single-layer 
solution  tested  earlier,  the  formulas  are  considerably  more  complicated,  and  the  eigenvalues  are  very  hard  to 
find.  The  computer  code  used  to  get  this  solution  is  therefore  very  complex;  it  will  be  discussed  in  more 
detail  in  Chapter  5.  For  now,  discussion  will  be  limited  to  a  comparison  with  a  simpler  solution. 

A  test  of  the  initial  condition  would  be  trivial  at  this  point;  inspection  of  Eq.  (4.75)  shows  that  every 
term  in  the  series  will  be  zero  for  f  =  0 ,  because  of  the  removal  of  the  term  F[X„) ,  in  accordance  with  the 
initial  condition.  The  temperature  will  be  zero  for  every  value  of  y,  a  test  would  yield  correct  values,  but 
this  would  not  say  much  about  the  validity  of  the  derivation.  However,  as  t  approaches  infinity,  the  steady- 
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state  temperatures  should  serve  as  a  check  of  the  eigenvalues,  the  eigenfunction  constants  given  in  Eqs. 
(4.71),  and  the  orthogonality  constant  given  by  Eq.  (4.82). 

Removing  the  generation  term  and  time  dependency  from  the  governing  equations  (4.59)  will  give  a 
simple,  linear,  steady-state  temperature  distribution.  This  new  problem  is  stated  as: 


(4.85a) 


(4.85a) 


with  the  boundary  conditions 


57,(0) 


+  //,r,(o)  =  //,7„ 


(4.86a) 


^5r,(L,)_57,(L,) 

dy  dy 


(4.86b) 


^^  =  W2[7,(L,)-72(7,)] 


(4.86c) 


dT^) 


+  Hyh{L,)  =  H,T, 


(4.86d) 


Note  that  because  this  problem  is  not  time-dependent,  neither  specific  heat  nor  the  thermal  diffusivity 
appears  in  the  equations.  The  governing  equations  (4.85)  are  integrated  twice,  resulting  in: 


7,  -  C,  y  -I-  D, 


(4.87a) 


Tj  =C2y  +  D2 


(4.87b) 


These  general  solutions  are  substituted  into  the  boundary  conditions  (4.86),  and  the  coefficients  C  and  D  are 
solved  for  through  some  simple,  if  tedious,  algebra. 


_ T]zL _ ^ 


(4.88a) 
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=?:+- 


Jr-L 


1  +  //, 


c,=- 


~+Ly-KL,+KL,+  — 
Ti-T^ 


KH, 


^KH, 


+  KL^-L^  +  L^+-^ 


{4.88b) 


(4.88c) 


D,=T,- 


a:l,+ 


I  K 

+  L^-KL,+KL,  +  — 
■ 


(4.88d) 


This  solution  can  be  used  to  check  the  two-layer  solution  by  substituting  a  sufficiently  large  value  of 
time  in  the  time-dependent  solution,  thus  allowing  it  to  approach  steady  state.  The  generation  term  is  set  to 
zero;  the  term  containing  it  in  Eq.  (4.76a)  becomes  zero  and  is  ignored. 

Table  4.4  shows  the  variables  used  to  calculate  the  temperature  distributions  for  the  two  models.  The 
physical  properties  for  layer  one  are  those  of  PMMA;  all  other  values  were  chosen  to  make  the  match  of  the 
two  solutions  more  apparent. 


Table  4.4.  Constants  for  test  of  two-layer  solution  with  no  generation 


0.05  cm 

t 

1.0  hr 

0.198  * 

0.10  cm 

h, 

4x10"-^ 

m  K 

1-0* 

'h 

1000 

iirK 

10.0  K 

3xl0'’-iip 

3x10"-^ 

K 

T, 

2.0  K 

The  temperature  distributions  are  given  in  Fig.  4.3.  The  two  solutions  agree  well  for  most  interior 
points,  with  the  infinite  sum  in  the  time-dependent  solution  truncated  after  150  terms.  The  agreement  of  the 
time-dependent  solution  at  the  endpoints  is  quite  bad,  because  of  the  truncation  of  the  summation.  Each  of 
the  two  governing  equations  is  defined  at  the  boundary  point  y  =  L, .  The  discontinuity  in  the  temperature 
at  L|  does  not  signify  a  problem  with  the  solution;  there  is  some  resistance  to  conduction  at  the  interface 
(because  of  the  small  value  of  /ij ),  and  the  temperature  should  be  expected  to  experience  a  jump  there. 
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Fig.  4.3.  Temperature  distribution  for  steady  state  test  of  two-layer  solution 

Table  4.5  shows  the  temperatures  for  the  resist  and  substrate  solutions  at  the  boundary  points.  There 
IS  excellent  agreement  at  the  boundary  between  the  resist  and  the  substrate,  but  the  time-dependent  solution 
IS  off  by  more  than  25%  at  the  substrate  surface.  This  should  not  detract  from  the  usefulness  of  the  solution. 

A  good  approximation  of  the  correct  temperature  may  be  inferred  from  a  graph  of  the  data  such  as  that  in 
Fig.  4.3. 


Table  4.5.  Boundary  point  temperatures  for  steady  state  test 


Surface: 

Temperature  Rise  (‘’C) 

Time  dependent 

Steady  state 

top 

8.520 

9.951 

boundary  ( r, ) 

5.009 

5.004 

boundary  ( ) 

3.080 

3.058 

bottom 

1.424 

2.066 

4.3.1.  Generation  Terms 

The  generation  function  for  this  solution  appears  in  the  coefficient  term  of  the  general  solution,  Eq. 
(4.76a).  Specifically,  the  term  we  are  now  interested  in  is 
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S  L  J  )dy 

7=1 

8j{y'^^')  given  by  Eqs.  (4.62).  Substituting,  we  get 


7-1 


where  Fj  is  L|  for  7  -  1 ,  and (Z^  +  for  7  =  2 ,  to  agree  with  the  definitions  in  Chapter  3. 
Substituting  the  eigenfunction  (4.67)  into  Eq.  (4.90)  and  rearranging  gives: 


7=1 


A.m  sin +  B.,„  cos^fL 
a,  /a,  ^ 


Integrating  this  expression  gives 


(4.89) 


(4.90) 


(4.91) 


7=1 


ju  sin-^^^ — ^cos-^’ 
a,.  Ja,  Ja 


A 


+  ^j.m 


•  J 


^n,y'  .  A.„V 

/ryCos-^+-^sin 

a. 


«y 


(4.92) 


Substituting  the  appropriate  limits  and  variables,  and  using  the  simplifications  given  in  Eq.  (4.69): 

^;;?^[A.™(/i.sinyL,-ycos7Z^  +  ye-'''"')  +  fi,„(/y,cosyL,  +  ysinyL,-/t,c-'''"')j+  (4.93) 

sin  77 -  T] COST)  L,)  (/y^ sin r]L,  -  rjcosriL, )J  + 

(^2  cosT)Lj  +  77sin77L3)-  e-''2'-2  (//^  cost)L,  +  7)sin  77L,)]} 

This  solution  is  substituted  into  Eq.  (4.76a).  The  derivation  is  completed  by  observing  that  with  an  initial 
condition  of  zero  temperature  rise  in  the  wafer,  the  functions  /^.(y)  are  zero,  making  the  function  F(A^) 
zero  also.  This  initial  condition  therefore  makes  Eq.  (4.74b)  equal  to  zero,  and  the  derivation  of  the 
solution  is  complete. 

4.3.2.  Confirmation  of  the  Generation  Terms 

The  two-layer  solution  with  generation  may  be  compared  to  the  single-layer  solution.  This 
comparison  will  help  to  verify  the  derivations  of  the  generation  terms  in  both  solutions.  With  no  generation 
in  the  substrate  and  very  high  conductivity  in  that  layer,  the  temperature  distribution  in  the  resist  is 
compared  with  a  cross  section  of  the  temperature  for  the  one-layer  problem.  To  check  that  the  generation  in 
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the  substrate  is  correct,  the  variables  for  the  two-layer  problem  are  simply  reversed;  the  temperature  profile 
should  be  reversed  also. 

The  constants  used  for  the  first  test  are  given  in  Table  4.6.  The  beam  width  for  the  two-dimensional 
solution  is  the  width  of  the  wafer,  and  the  velocity  is  almost  zero;  this  approximates  the  one-dimensional, 
two-layer  solution.  A  relatively  low  value  of  time  is  used,  to  minimize  the  effects  of  the  end  conditions  on 
the  two-dimensional  solution.  The  high  values  of  ,  A3 ,  and  a  2  make  the  substrate  in  the  two-layer 
solution  essentially  disappear;  that  is,  the  resistance  to  heat  transfer  is  so  small  that  the  boundary  condition 
at  the  interface  in  the  two-layer  solution  is  essentially  identical  to  the  zero-temperature  boundary  condition 
on  the  one-layer  problem. 


Table  4.6.  Constants  for  test  of  the  resist  generation  term 


One-layer  solution 

L  0.03  cm 

W  -1.082x10’-^ 

nr 

k  0.198  * 

a  1 00.0  cm 

-2310^ 

h  600 

in'K 

a  I.i8xl0‘''4 

w  95.0  cm 

vO.OOOl  IxlO-^f 

Two-layer  solution 

L,  0.03  cm 

a  6.0  cm 

t  4.0  s 

Lj  0.1  cm 

W,  2.164  X 10’ 

*  in 

k,  0.198  ^ 

h.  10 

1  m  K 

0 

d 

k,  IxlO’* 

fi,  2310^ 

a,  1.18x10'’ -if 

/X,  lOO.Oi 

ttj  lxl0“*4 

The  results  of  this  test  are  shown  in  Table  4.7;  the  two  solutions  agree  very  well,  with  small 
differences  due  to  convergence  and  computational  errors. 


Table  4.7.  Temperatures  (in  °C)  from  test  of  the  two-layer  solution 


y  Position  (cm) 

One-Layer 

Two-Layer 

Normal 

Reversed 

Surface 
Midpoint 
Interface  or  Bottom 

10.6403 

8.2623 

1.3e-14 

10.6426 

8.2627 

0.00016 

10.6420 

8.2631 

0.00024 
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To  test  the  substrate  generation  term,  the  constants  were  redefined  to  reverse  the  orientation  of  the 
wafer.  The  convection  coefficients  and  li^  are  simply  swapped,  as  are  the  conductivities  and  thermal 
diffusivities.  The  irradiation  and  absorption  coefficient  must  be  redefined  because  of  the  form  of  the 
generation  terms.  The  interface  conduction  coefficient  h2  remains  the  same.  Although  the  physical 
situation  has  been  reversed,  a  simple  change  of  variables  doesn’t  transform  the  solution  to  an  exact  reverse 
of  the  original.  For  instance,  the  generation  term  is  defined  for  an  energy  source  which  decays 
exponentially  with  depth.  The  constants  Wj  and  defined  to  reverse  the  decay  profile  for  the 

substrate,  but  the  new  values  are  not  simply  the  opposite  of  the  originals.  Because  of  this  nonsymmeU'y, 
truncation  and  roundoff  error  will  introduce  small  variations  in  the  temperatures,  as  shown  in  Table  4.7. 


CHAPTER  5 


ANALYSIS  OF  SOLUTIONS 

Solving  the  differential  equations  of  conduction  heat  transfer  is  only  the  first  step  in  getting  a 
temperature  distribution.  For  series  solutions  such  as  the  ones  presented  in  Chapter  4,  evaluation  of  the 
series  can  be  more  difficult  than  the  mathematics  involved  in  obtaining  it. 

In  this  chapter,  some  of  the  characteristics  of  series  solutions  will  be  discussed,  and  the  impact  of  these 
characteristics  on  the  evaluation  of  solutions  will  be  examined.  The  computer  programs  used  to  evaluate 
the  solutions  are  discussed,  along  with  an  analysis  of  the  programs’  advantages  and  limitations.  Some 
temperature  distributions  produced  by  the  computer  programs  are  presented  for  a  variety  of  test  cases. 

Experimental  tests  ot  the  temperature  rises  were  performed  at  the  Center  for  Advanced 
Microstructures  and  Devices  in  Baton  Rouge.  The  chapter  concludes  with  an  explanation  of  the 
experimental  procedure  and  the  results  of  several  exposure  tests. 

5.1.  Series  Solutions 

The  most  basic  problem  encountered  in  evaluation  of  series  solutions  is  that  the  infinite  series  must  be 
truncated  at  some  point.  This  error  may  be  minimized  by  calculating  a  large  number  of  terms  in  the  series. 
For  many  series,  though,  this  is  not  necessary;  only  a  few  terms  may  be  adequate.  In  general,  the  magnitude 
of  the  terms  in  the  series  decreases  as  the  index  (here  denoted  as  m  or  n)  increases,  and  this  decrease  may  be 
approximately  linear,  or  may  be  some  other  function  of  the  index.  For  the  solutions  presented  in  Chapter  4, 
the  eigenvalues  are  complex  functions  of  the  indices,  and  the  terms  themselves  are  very  complicated.  Some 
experimentation  is  therefore  necessary  to  find  the  number  of  terms  needed  for  an  accurate  solution. 

The  shape  and  complexity  of  the  function  being  represented  can  give  an  indication  of  what  the 
convergence  of  its  series  will  be.  For  instance,  the  surface  temperature  for  the  one-layer,  two-dimensional 
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solution  has  a  relatively  complicated  shape,  which  bears  very  little  resemblance  to  the  sine  function  which  is 
used  to  represent  it  in  the  series.  The  series  therefore  converges  slowly  (see  Fig.  5.1).  Summing  70  terms 
gives  a  very  good  solution,  and  ten  terms  is  a  close  approximation,  with  sharp  edges  rounded.  Four  terms 
gives  a  very  bad  approximation.  Compare  this  to  Fig.  5.2,  which  shows  the  convergence  of  the  same 
solution  in  the  y  direction.  Because  of  the  simplicity  of  the  temperature  profile,  one  term  gives  a  reasonable 
approximation,  and  any  summation  of  more  than  four  terms  is  essentially  the  same.  The  number  of  terms 
summed  for  each  of  the  profiles  in  these  two  graphs  is  indicated  by  n  and  m,  respectively. 


Another  way  of  characterizing  the  convergence  is  by  examining  the  magnitude  of  the  individual  terms 
in  the  series.  Each  term  consists  of  sine  and/or  cosine  waves,  and  the  maximum  value  of  each  term  may  be 
used  as  a  measure  of  the  importance  of  that  term.  The  error  in  the  final  solution  depends  on  where  the 
series  is  truncated,  and  on  the  magnitudes  of  those  terms  which  are  not  added  to  the  solution.  The  poor 
convergence  of  the  two-layer  solution  is  illustrated  by  the  scatter  graph  shown  in  Fig.  5.3. 
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y  (urn) 

Fig.  5.2.  Convergence  of  one-layer  solution  in  y  direction 
at  the  center  of  the  resist 
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Fig  5.3.  Convergence  of  one  and  two  layer  solutions 


The  values  for  the  one-layer  solution  in  Fig.  5.3  represent  the  maximum  value  of  the  term  for  any 
particular  y-direction  eigenvalue;  no  particular  .x-direction  eigenvalue  is  used.  The  convergence  of  the  one- 
layer  solution  is  very  uniform  -  the  magnitude  of  the  terms  decreases  quickly  and  predictably.  By  the  tenth 
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term,  there  is  a  decrease  of  three  orders  of  magnitude.  This  explains  the  accuracy  of  the  limited  summations 
shown  m  Fig.  5.2.  The  resist  and  substrate  summations  for  the  two-layer  solution  do  not  converge  well  at 
all.  By  the  one-hundredth  term,  there  is  a  decrease  of  only  two  orders  of  magnitude  in  the  terms.  For  these 
summations,  the  magnitude  does  not  decrease  in  any  predictable  fashion;  the  data  points  in  Fig.  5.3  are 
spread  over  a  range  of  three  or  four  orders  of  magnitude  for  successive  terms. 

5.2.  Computer  Programs 

The  solutions  presented  herein  were  evaluated  using  programs  written  by  the  author.  Although  there 
are  commercially  available  programs  that  are  capable  of  evaluating  these  series,  there  are  two  reasons  why 
these  were  not  used.  First,  because  of  their  generality  and  user  interfacing,  these  programs  are  quite  slow, 
and  use  a  lot  of  memory.  Mathcad  5.0  was  able  to  evaluate  some  of  the  preliminary  solutions,  but  took 
several  minutes  to  calculate  abbreviated  sums.  The  number  of  terms  evaluated  was  severely  limited  by  the 
memory  of  the  machine. 

Time  delays  are  simply  a  nuisance;  the  primary  reason  that  programs  were  written  from  scratch  was  to 
insure  precise  knowledge  and  control  of  the  numerical  methods  used  in  the  evaluation  of  the  solution.  This 
was  particularly  important  in  the  finding  of  eigenvalues,  because  of  the  transcendental  equations  involved. 
5.2.1.  Finding  Eigenvalues 

Formulas  such  as  Eq.  (4.13)  are  easily  evaluated,  but  solving  transcendental  equations  such  as  Eq. 
(4.19)  or  Eq.  (4.72)  can  be  quite  complicated,  and  it  is  extremely  important  that  all  of  the  desired 
eigenvalues  are  found.  The  first  one  is  often  the  most  important,  and  the  hardest  to  find. 

For  the  one  layer  problem,  the  eigenvalues  are  the  roots  of  the  equation 

~^„co\.k„L=H  (5.1) 

Specifically,  the  program  finds  the  zeros  of  the  function 

/(Aj  =  //  +  A„cotA„L  (5.2) 

The  routine  used  to  evaluate  Eq.  (5.1)  is  relatively  simple,  and  very  reliable,  because  the  eigenvalues  lie 
within  intervals  of  equal  length,  bounded  by  singularities.  The  positions  of  these  singularities  correspond  to 
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the  period  of  the  cotangent  function  in  Eq.  (5.2).  Figure  5.4  shows  a  graph  of  this  function.  The  vertical 
lines  show  the  positions  of  the  singularities,  at  which  the  function  approaches  positive  or  negative  infinity. 


Fig.  5.4.  Graph  of  the  eigenvalue  function  for  one-layer  problem 

Because  of  the  nature  of  the  cotangent  function,  an  intersection  point  must  lie  in  each  interval.  The 
positions  of  the  singularities  are  calculated  from  the  material  properties  and  dimensions  (which  determine 
the  argument  in  the  cotangent  function),  and  Bidder’s  method  is  used  to  find  the  roots  (Press  et  al.  1992). 
This  is  a  simple  and  efficient  routine,  which  is  perfectly  suited  to  this  purpose;  given  two  points  bracketing 
a  root,  the  method  will  find  that  root,  without  any  possibility  of  jumping  out  of  the  interval  to  settle  on  a 
different  one.  It  also  has  the  advantage  of  converging  very  quickly,  gaining  as  much  as  two  or  three 
significant  digits  of  accuracy  with  each  iteration. 

Evaluation  of  the  eigenvalues  for  the  two-layer  problem  proved  to  be  much  more  difficult.  These 
eigenvalues  are  given  by  the  zeros  of  the  determinant: 


r 

-H,  0 

0 

.  ysinil, 

cos>l,-i-  ^  '  -sinrjL, 

"2 

-costjL, 

Kyco%yL^ 

-Ky  Sin  yL^  -7]cost?L, 

r;sin  rjL, 

0 

0  "'“"‘’.SMI, 

^3 

7]sinnL, 

+cosr7L2 

With  all  the  dependencies  and  different  functions  in  this  transcendental  equation,  no  simple  means  of 
finding  the  roots  could  be  devised.  The  appearance  of  the  determinant  may  vary  greatly,  depending  on 
constants  and  dimensions  used,  but  a  typical  graph  is  given  in  Fig.  5.5. 


Fig.  5.5.  Determinant  of  the  eigenvalue  matrix  for  two-layer  problem 


This  function  is  highly  irregular,  and  the  roots  do  not  lie  on  any  obvious  intervals,  making  them  much 
harder  to  find.  The  position  and  spacing  of  the  roots  is  dependent  on  a  variety  of  constants,  and  will  be 
different  for  any  two  cases.  This  requires  an  algorithm  that  is  relatively  flexible,  yet  reliable,  allowing  it  to 
be  used  for  a  variety  of  cases  without  having  to  modify  the  search  routine  and  recompile  the  code. 

Fortunately,  the  determinant  is  finite,  at  least  over  the  range  of  interest.  This  means  that  a  simple  step 
routine  can  be  use  to  find  the  roots.  The  routine  works  by  evaluating  the  function  at  the  beginning  and  end 
of  an  interval  (determined  by  the  step  size),  and  checking  the  sign  of  the  values  to  see  if  the  function  has 
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crossed  the  axis.  A  sign  change  over  an  interval  always  means  that  a  root  lies  on  that  interval  (with  the 
single-layer  eigenvalue  function,  this  may  have  signified  a  singularity).  Once  a  root  has  been  bracketed  in 
this  manner,  Ridder’s  method  (Press  et  al.  1992)  can  be  used  to  pinpoint  it.  The  difficulty  with  this 
algorithm  is  that  if  the  step  is  too  large,  it  might  miss  roots  altogether.  If  the  function  has  crossed  the  axis 
twice  within  the  interval,  the  signs  of  the  endpoint  values  will  be  the  same,  and  the  algorithm  will  not 
recognize  that  a  root  exists  on  the  interval. 

The  routine  used  in  subroutine  LAMDA  of  program  THESIS3.CPP  in  the  appendix  uses  a  variable 
step  size,  which  was  refined  during  evaluation  of  various  solutions.  Because  the  first  root  is  generally  very 
close  to  the  origin,  the  step  size  starts  off  small.  After  the  first  root  is  found,  the  step  is  increased,  based  on 
the  distance  between  previous  roots.  This  routine  has  worked  without  error  for  a  variety  of  cases;  however, 
if  the  program  gives  unexpected  results,  the  eigenvalues  should  be  checked,  as  the  method  is  by  no  means 

foolproof.  This  may  be  easily  done  by  graphing  the  eigenvalue  function  (as  shown  in  Fig.  5.5)  and 
estimating  the  eigenvalues. 

5.2.2.  Series  Summation 

Finding  the  temperature  at  even  one  point  in  space  requires  a  large  number  of  calculations.  For  the 
two-dimensional  solution,  taking  only  thirty  terms  in  each  summation  requires  that  many  of  the  formulas  in 
the  terms  be  calculated  nine  hundred  times.  For  the  multiple-pass  solution,  finding  the  temperature  at  one 
point  at  a  specific  time  requires  the  evaluation  of  the  two-dimensional  solution  for  several  passes. 
Obtaining  the  temperature  at  several  points  as  a  function  of  time  requires  that  many  terms  be  calculated 
thousands  of  times,  with  the  terms  themselves  containing  several  calculations.  This  was  unacceptable, 
especially  for  testing  and  debugging  purposes,  when  many  cases  were  being  studied. 

One  way  to  achieve  significant  time  savings  is  to  recognize  that  the  position  variables  x  and  y  appear 
only  m  the  arguments  of  the  sine  and  cosine  terms  of  the  eigenfunctions.  All  other  calculations  in  the  terms 
are  the  same  for  every  point,  and  need  only  be  performed  once.  The  term  is  then  completed  for  each  point 
by  multiplying  by  the  sine  or  cosine  term  containing  x  or  y,  and  adding  this  to  the  running  total  of  the 
summation  for  the  respective  point.  Although  the  temperatures  must  be  stored  in  an  array  because  they  are 
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solved  simultaneously,  the  total  number  of  calculations  is  greatly  reduced.  This  in  turn  reduces  the  time 
needed  to  get  a  temperature  distribution,  even  for  a  one-dimensional  solution  with  a  small  number  of  points. 
For  a  two-dimensional,  multiple-pass  solution,  the  time  savings  is  even  larger. 

Another  way  to  cut  down  on  the  number  of  calculations  is  to  take  as  many  of  them  as  possible  out  of 
the  loops  for  the  summations.  For  instance,  the  eigenvalues  arc  the  first  thing  to  be  removed  from  the  loops; 
their  calculation  is  very  time  consuming,  and  they  are  saved  into  arrays  before  the  summation  loops  are 
begun.  This  means  that  the  eigenvalues  need  only  be  calculated  once,  and  storing  them  in  arrays  does  not 
use  much  of  the  computer’s  memory.  Other  parts  of  the  terms  which  depend  only  on  one  set  of  eigenvalues 
may  be  separated  from  one  of  the  summations.  For  instance,  the  calculations 

4a  M'/i ,  ,  ,  , 

~7Z  ~A:  ^m)  (5.4) 

can  be  removed  from  the  term 


Cosfi„{tv  +  w)~cosl}^tv-e  1)] 

sin  (rv  +  iv) -sin  Pjv - e~ sin P„\ 


(5.5) 


Pn'' 


which  is  the  coefficient  of  the  two  dimensional  solution  given  by  Eq.  (4.35).  The  group  of  terms  in  Eq. 


(5.4)  can  then  be  taken  out  of  the  loop  for  )3„ ,  so  the  calculations  are  done  only  once  for  each  in  term, 

instead  of  being  done  in  each  of  the  fifty  or  more  n  terms  needed  to  get  an  accurate  value  of  temperature  for 
this  solution. 


5.3.  Test  Cases 

Some  observations  may  be  made  about  the  expected  temperature  distributions  and  resulting 
deformations.  Strictly  speaking,  there  will  never  be  a  steady-state  temperature  distribution  in  the  wafer, 
because  the  X-ray  beam  is  always  moving  with  respect  to  the  wafer.  However,  after  a  number  of  passes 
have  heated  the  wafer,  the  heat  generation  from  the  beam  will  equal  the  total  dissipation  through  convection 
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to  the  surroundings  and  through  conduction  into  the  substrate.  The  time  previous  to  this  condition  might  be 
called  a  warm-up  time.  Deformations  in  the  resist  will  presumably  be  relatively  large  when  going  from  time 
zero  to  the  warm-up  time,  then  drop  to  smaller  values  as  the  temperature  fluctuates  with  each  pass.  Being 
able  to  establish  the  length  of  this  warm-up  time  and  the  magnitude  of  the  temperature  differences  involved 
will  allow  the  researcher  to  determine  whether  exposure  during  these  first  few  passes  will  have  a  detrimental 
effect  on  the  final  structure. 

5.3.1.  Resist  Deformations 

Thermal  stresses  in  the  wafer  along  the  boundary  between  the  resist  and  the  substrate  may  weaken  the 
bond  between  the  two  materials.  Bonding  problems  are  often  observed,  so  one  should  expect  to  find 
significant  stresses  along  the  interface.  A  two-layer,  two-dimensional  solution  would  have  provided  the 
best  estimate  of  the  temperature  distributions  in  the  resist  and  substrate,  but  the  simpler  solutions  derived  in 
Chapter  4  can  give  some  insight  into  the  temperatures  and  resulting  stresses. 

Most  substrate  materials  have  high  thermal  conductivity,  so  the  temperature  in  the  substrate  should  be 
expected  to  be  almost  uniform,  and  thermal  gradients  in  the  substrate  itself  should  not  be  high  enough  to 
produce  any  deformation  of  the  resist.  The  two-layer,  one-dimensional  model  may  be  used  to  confirm  this 
theory.  Table  5.1  shows  the  constants  used  for  this  calculation.  The  intensity  of  the  X-ray  beam 
approximates  a  limiting  case;  it  is  as  high  as  one  might  expect  to  see  in  an  exposure  station  at  CAMD.  The 
heat  transfer  coefficients  are  also  chosen  for  a  limiting  case;  the  convection  on  the  resist  surface  is  zero, 
while  the  convection  on  the  back  of  the  substrate  is  high;  this  will  force  more  of  the  heat  to  flow  through  the 
substrate,  causing  a  higher  temperature  gradient  there. 

Figure  5.6  shows  the  temperature  distribution  in  the  wafer  after  20  seconds.  The  temperature  gradient 
across  the  substrate  is  only  2  °C,  while  the  gradient  across  the  resist  is  140  °C.  No  thermal  stresses  should 
be  expected  in  the  substrate  due  to  temperatures  within  the  substrate  itself.  This  is  not  to  say  that  there 
would  be  no  stresses  at  all;  but  they  would  be  caused  by  different  expansion  of  the  resist  and  substrate,  and 
the  temperature  of  the  substrate  may  safely  be  treated  as  uniform. 


Table  5.1 .  Constants  for  test  of  the  substrate  temperature  gradient 
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Fig.  5.6.  Temperatures  from  test  of  the  substrate  temperature  gradient 

Note  that  the  high  resist  temperatures  shown  in  Fig.  5.6  should  not  be  seen  in  practice.  After  20 
seconds,  the  one-dimensional  model  is  approaching  steady-state,  and  the  heat  transfer  coefficient  on  the 
resist  surface  is  set  to  zero,  but  should  be  about  600 

nr-K  * 

5.3.2.  Multiple  Pass  Temperatures 

The  scanning  of  the  wafer  through  the  X-ray  beam  will  cause  some  thermal  cycling  in  the  resist.  This 
will  probably  compound  any  stress-related  problems  in  the  exposure  and  developing  processes.  The  cycling 
may  cause  some  fatigue  in  the  bond  between  resist  and  substrate,  and  may  compound  the  problem  of 
misalignment  between  the  mask  and  the  wafer. 
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Two  types  of  misalignment  problems  may  be  experienced  during  the  manufacturing  process.  In  the 
first,  the  resist  is  deformed  by  some  steady-state  temperature,  remaining  in  one  position  for  most  of  the 
exposure.  When  the  wafer  is  cooled  after  exposure,  the  structures  imprinted  on  it  can  shift.  Edge  definition 
should  be  good,  though  the  structures  may  be  slightly  warped;  for  instance,  walls  may  not  be  straight.  The 
second  misalignment  problem  is  related  to  the  thermal  cycling.  If  the  resist  is  constantly  moving  during 
exposure,  some  areas  may  not  be  fully  exposed,  and  therefore  will  not  develop  properly  during  processing. 
Because  of  the  small  feature  dimensions  which  may  be  made  by  the  LIGA  process,  even  very  small 
temperature  gradients  may  cause  misalignment  problems.  An  analysis  of  the  temperature  fluctuations  in  the 
resist  as  a  function  of  time  may  give  some  insight  into  both  types  of  misalignment  problems. 


Table  5.2.  Constants  for  calculation  of  multiple 
pass  temperatures 
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Figure  5.7  shows  temperatures  with  respect  to  time  for  six  locations  in  the  resist,  calculated  with  the 
values  given  in  Table  5.2.  The  thick  gray  lines  trace  the  temperature  for  locations  at  the  center  of  the  resist; 
the  X-ray  beam  passes  over  these  points  at  regular  intervals.  The  thin  black  lines  trace  temperatures  near 
the  end  of  the  beam’s  travel;  when  the  beam  reverses  direction  at  that  end  of  the  resist,  it  passes  over  these 
points  twice  in  quick  succession,  resulting  in  a  higher  maximum  temperature  than  is  observed  at  the  center 
of  the  resist.  The  end  of  the  resist  then  has  almost  the  entire  period  of  oscillation  to  cool  down;  it  also 
reaches  a  lower  temperature  than  is  seen  at  the  center  of  the  resist.  Obviously,  the  thermal  cycling  of  the 
resist  can  vary  greatly  from  one  physical  location  to  another. 
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The  temperatures  and  the  magnitude  of  their  changes  are  smaller  near  the  interface  between  the  resist 
and  substrate,  compared  to  the  temperatures  at  the  surface.  This  should  result  in  smaller  local  stresses  near 
the  interface.  Near  the  surface,  tliere  will  be  higher  stresses  and  more  deformation;  it  is  here  that 
incomplete  development  will  be  a  problem.  Deformation  can  move  a  part  of  the  resist  with  respect  to  the 
mask;  this  misalignment  will  cause  parts  of  the  PMMA  near  the  edge  of  exposed  areas  to  move  in  or  out  of 
the  beam,  resulting  in  an  incorrect  dosage  and  unpredictable  structure  shapes. 

5.4.  Experimental  Results 

In  order  to  get  an  idea  of  the  actual  temperature  rises  in  thick  resists,  some  experiments  were 
performed  at  CAMD  in  Baton  Rouge.  Wafers  were  constructed  with  thermocouples  on  the  surface  and  at 
the  interface  between  resist  and  substrate.  The  temperatures  in  these  locations  were  then  measured  during 
exposure,  using  a  computer  and  data  acquisition  system. 
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Fig.  5.8.  Wafer  used  in  experimental  test  of  heat  transfer 


The  resist  used  in  these  experiments  was  0.5  mm  thick  PMMA,  and  the  substrate  was  a  0.5  mm  silicon 
wafer.  Dimensions  are  shown  in  Fig.  5.8.  The  inner  circle  in  the  figure,  labeled  “exposed  area,”  matches 
the  location  of  the  hole  in  the  mask  through  which  parts  of  the  X-ray  beam  (rectangle)  can  pass.  The 
thermocouples  were  J  type,  made  with  75  pm  iron  and  constantan  wires.  Though  very  small,  the  junction 
(where  the  temperature  is  measured)  was  quite  large  when  compared  to  the  500  pm  resist  and  substrate;  for 
this  reason,  the  junctions  were  hammered  flat  to  minimize  their  disruption  of  the  heat  transfer.  The  silicon 
wafers  used  in  the  making  of  microstructures  are  typically  coated  with  a  metal,  used  as  a  base  for 
electroplating  later  in  the  process.  Uncoated  silicon  was  used  in  this  experiment  to  avoid  shorting  the 
thermocouples  which  are  in  contact  with  the  substrate.  Because  of  the  high  conductivity  of  nickel,  and  the 
thinness  of  the  metal  layers  generally  used  in  this  application  (about  1  nm),  the  absence  of  this  layer  should 
not  significantly  change  the  heat  transfer  characteristics  of  the  wafer.  Silicon  is  a  semiconductor,  but 
preliminary  tests  showed  that  contact  with  the  silicon  wafer  had  no  measurable  effect  on  the  accuracy  of  the 
thermocouples.  The  small  size  of  the  junctions  produced  very  fast  response  times  in  the  thermocouples. 
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Two  wafers  were  made  for  this  experiment;  on  the  first,  the  bond  between  the  PMMA  and  the  silicon 
was  not  very  good,  with  small  air  pockets  under  much  of  the  PMMA.  This  might  produce  a  high  thermal 
resistance  at  that  junction,  causing  higher  temperatures  in  the  resist.  The  second  wafer  had  a  very  good 
resist-substrate  bond. 

The  first  irradiation  test  was  designed  around  a  worst-case  scenario.  A  wire  mesh  mask  was  used; 
about  80  percent  of  the  radiation  incident  on  the  mask  passes  through  to  the  wafer.  The  positions  of  the 
thermocouples  are  shown  in  Fig.  5.8.  The  synchrotron  was  running  at  1.5  GeV,  and  between  125  and  100 
mA,  and  the  X-ray  beam  was  not  filtered  before  the  irradiation  chamber.  The  scan  length  was  3.8  cm,  with 
a  scan  velocity  of  0.635  cm/s,  and  the  pressure  in  the  chamber  was  10  torr.  A  500  micron  layer  of  PMMA 
cannot  be  successfully  exposed  under  these  conditions.  After  30  minutes,  the  resist  was  overexposed,  and 
the  surface  of  the  PMMA  appeared  puffy  and  discolored,  as  if  bubbled  from  heat.  Figure  5.9  shows  the 
temperature  ranges  near  the  beginning  and  the  end  of  the  run.  The  maximum  temperature  measured  in  the 
resist  was  50°C,  at  thermocouple  3. 

During  the  first  test,  thermocouple  3  was  inadvertently  placed  Just  inside  the  exposed  area  of  the  resist, 
so  the  temperatures  from  that  location  are  comparable  to  the  ones  measured  at  thermocouple  4.  This 
supports  the  assumption  made  elsewhere  that  variations  in  the  x  direction  (parallel  to  the  surface)  are  small 
compared  to  variations  in  the  y  direction  (normal  to  the  surface). 

Because  the  measured  temperatures  were  so  low  (less  than  half  of  the  melting  point  of  PMMA),  the 
bubbling  of  the  PMMA  must  be  assumed  to  be  caused  by  a  chemical  phenomenon,  and  not  by  high 
temperatures  in  the  resist.  The  chemical  reaction  caused  by  irradiation  of  PMMA  with  X-rays  may  produce 
gasses  as  a  by-product.  It  appears  that  the  high  dose  is  preventing  these  gasses  from  escaping  quickly 
enough,  causing  them  to  swell  the  PMMA. 
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Fig.  5.9.  Temperature  ranges  for  the  first  exposure  test 


Thernocouple 


Fig.  5.10.  Temperature  ranges  for  the  second  exposure  test 


For  the  second  and  third  exposures  of  the  temperature  test,  the  better  wafer  was  used.  A  14  |xm 
aluminum  filter  was  placed  in  the  beam  line,  and  at  the  beginning  of  the  second  exposure  the  current  in  the 
ring  had  fallen  to  about  80  mA,  at  1.5  GeV.  All  other  parameters  were  unchanged.  These  conditions  are 
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similar  to  those  of  a  typical  exposure.  Figure  5.10  shows  the  temperature  ranges  for  this  test  run.  The 
temperatures  are  quite  low,  with  the  fluctuations  at  the  center  limited  to  about  33°C. 

Before  the  third  run,  the  synchrotron  was  reinjected,  increasing  the  current  to  about  130  mA.  The 
wafer  was  also  turned  90  degrees  and  shifted  such  that  the  thermocouples  were  near  the  top  and  bottom  of 
the  area  exposed  to  the  X-ray  beam;  that  is,  they  were  near  the  ends  of  the  beam’s  travel. 


Thernocouple 

Fig.  5.11.  Temperature  ranges  for  the  third  exposure  test 


Figure  5.1 1  shows  the  temperature  ranges  for  the  third  run.  The  increased  current  in  the  ring  probably 
caused  most  of  the  difference  between  these  temperatures  and  those  obtained  in  the  previous  exposure.  It  is 
worth  noting,  however,  that  the  temperature  range  in  exposure  three  is  a  much  larger  percentage  of  the 
maximum  temperature;  the  ends  of  the  exposed  area  have  almost  the  whole  period  of  the  beam  oscillation  to 
cool  down  between  passes,  while  the  center  of  the  wafer  has  only  half  that  amount  of  cool-down  time.  This 
is  the  same  phenomenon  that  was  illustrated  in  Fig.  5.7. 


5.5.  Conclusions 

Though  the  analytical  solutions  derived  in  Chapter  4  have  some  inherent  limitations,  such  as 
truncation  of  the  series  and  the  difficulty  of  finding  eigenvalues,  these  were  successfully  overcome  by 
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designing  computer  programs  specific  to  this  application.  The  programs  appeared  to  be  completely 
successful  in  finding  all  the  eigenvalues,  both  for  the  tests  presented  in  this  chapter  and  for  many  others  not 
included  here.  For  the  one-layer  solution,  truncation  does  not  appear  to  be  a  problem,  even  when  a 
moderate  number  of  terms  are  used.  However,  the  two-layer  solution  requires  a  large  number  of  terms  to 
insure  accuracy,  especially  near  the  surfaces  of  the  resist  and  substrate.  The  maximum  number  of  terms 

handled  by  program  THESIS4  (found  in  the  appendix)  is  300;  this  is  the  number  that  was  used  for  all  of  the 
tests  presented  here. 

The  analytical  and  experimental  tests  both  indicate  that  temperature  rises  in  the  resist  will  be  relatively 
small;  at  most,  one  might  expect  to  see  temperature  fluctuations  of  20°C  when  using  scan  speeds  of  at  least 
0.635  cm/s.  Most  masks  allow  only  a  small  percentage  of  the  X-ray  radiation  to  pass  through  to  the  resist, 
and  the  synchrotron  is  rarely  running  at  its  maximum  energy.  For  most  combinations  of  mask  and  wafer, 
the  temperature  rises  should  be  a  few  degrees,  at  most;  the  resist  deformations  produced  by  fluctuations  of 
this  magnitude  should  be  very  small. 


CHAPTER  6 


CONCLUSIONS  AND  RECOMMENDATIONS 

A  number  of  conclusions  may  be  drawn  from  the  analytical  and  experimental  work  presented  here,  not 
only  about  temperature  distributions  in  irradiated  wafers,  but  also  about  the  best  methods  of  finding  these 
temperatures,  and  the  usefulness  of  these  methods  in  the  study  of  deformations. 

The  analytical  and  experimental  tests  both  indicate  that  temperature  rises  in  the  resist  will  be  relatively 
small;  at  most,  one  might  expect  to  see  temperature  fluctuations  of  20°C  when  using  scan  speeds  of  at  least 
0.635  cm/s.  Faster  scan  speeds  will  produce  lower  maximum  temperatures  and  smaller  temperature 
fluctuations.  For  most  exposure  situations,  the  temperature  rises  should  be  a  few  degrees,  at  most.  The 
resist  deformations  resulting  from  such  small  temperature  fluctuations  should  be  minimal. 

6.1.  Analytical  Models 

Two  principle  analytical  solutions  were  developed  in  Chapter  4,  and  the  computer  code  necessary  to 
evaluate  them  was  presented  in  Chapter  5.  The  accuracy  of  the  temperatures  generated  by  these  programs 
with  respect  to  the  governing  equations  and  boundary  conditions  has  been  established.  However,  this  does 
not  guarantee  a  good  approximation  of  the  actual  temperature  distribution  during  exposure. 

A  number  of  simplifications  were  necessary  to  construct  problems  which  could  be  solved  analytically. 
It  was  discovered  that  multiple-layer  solutions  are  limited  to  one  dimension  except  in  special  cases  (see 
Chapter  3);  this  prevented  the  realization  of  a  two-dimensional,  two-layer  analytical  solution.  For  the  one- 
layer  problem,  the  rectangular  shape  of  the  beam  necessitated  a  Cartesian  coordinate  system.  A  zero 
temperature  boundary  condition  was  assumed  for  the  lower  surface,  and  a  convection  condition  was  used 
for  the  upper  surface,  which  is  in  fact  bordered  by  a  very  thin  layer  of  gas  and  a  mask;  this  mask  has  a 
temperature  distribution  of  its  own.  It  was  also  not  possible  to  take  into  account  the  shape  or  thermal  mass 
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of  the  carrier  which  holds  the  wafer  during  exposure.  The  utility  of  both  of  the  analytical  models  is  severely 
limited  by  the  large  number  of  simplifications. 

The  solutions  to  the  analytical  models  must  be  evaluated  computationally,  because  of  their  complexity. 
The  programs  used  for  evaluation  of  these  solutions  were  very  difficult  to  develop. 

Another  limitation  of  the  analytical  models  is  that  they  cannot  be  used  for  analysis  of  local 
deformations  in  the  resist.  The  generation  terms  can  only  model  the  part  of  the  X-ray  beam  that  passes 
through  the  mask  as  a  rectangle  of  energy;  heat  generation  in  the  shape  of  individual  microstructures  would 
require  new  and  far  more  complicated  generation  terms.  Bulk  deformations  will  undoubtedly  cause 
alignment  problems.  However,  with  certain  masks,  large  scale  temperatures  predicted  by  the  analytical 
models  might  be  insignificant,  while  local  temperature  changes  in  the  vicinity  of  small  irradiated  areas 
might  be  very  important.  This  phenomenon  cannot  be  adequately  handled  with  the  analytical  solutions 
presented  here,  and  construction  of  analytical  solutions  which  could  produce  data  of  this  type  would  be 
subject  to  all  of  the  simplification  problems  described  above. 

6.2.  Experimental  Analysis 

The  experiments  described  in  Chapter  5  produced  some  interesting  and  potentially  valuable  data, 
though  the  design  was  relatively  simple.  Its  primary  flaw  is  the  presence  of  the  thermocouples.  Metals 
absorb  much  more  radiation  than  PMMA;  if  the  size  of  a  thermocouple  junction  is  too  large,  it  might  absorb 
enough  radiation  to  raise  its  temperature  significantly  above  the  temperature  of  the  surrounding  PMMA.  A 
large  junction  may  also  interrupt  the  normal  flow  of  heat  in  the  wafer. 

With  only  four  thermocouples,  too  little  data  was  gathered  from  these  experiments  to  perform  a 
detailed  analysis  of  deformations  in  the  resist.  Though  other  methods  of  measuring  the  temperatures  may 
be  used,  it  would  most  likely  be  very  difficult  get  measurements  with  high  enough  spacial  resolution  for  an 
analysis  of  local  deformations  and  the  effect  of  those  deformations  on  the  shape  of  the  finished 
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6.3.  Recommendations 

The  usefulness  of  the  analytical  solutions  presented  herein  is  severely  limited  by  the  simplifications 
necessary  to  construct  solvable  problems.  In  theory,  these  simplifications  need  not  be  made  when  using  a 
finite  element  or  finite  difference  model.  Such  programs  can  more  easily  handle  complexities  such  as 
multiple  layers,  different  materials,  convection  heat  transfer,  and  unusual  generation  profiles.  If  one  is 
concerned  with  the  accuracy  of  a  numerical  solution,  keep  in  mind  the  complexity  of  the  programs  required 
for  the  evaluation  of  the  analytical  solutions  presented  here,  and  the  inaccuracies  produced  by  truncation. 
From  a  heat  transfer  standpoint,  there  is  nothing  particularly  unusual  about  the  wafer,  the  structure  that 
supports  it,  or  the  environment  in  which  it  is  irradiated,  which  should  be  challenging  to  model  numerically. 

A  numerical  model  of  the  resist  would  not  only  allow  a  more  reliable  study  of  the  large-scale 
temperature  changes  in  the  resist,  but  would  enable  the  researcher  to  study  more  complex  situations,  such  as 
the  temperature  field  and  deformations  around  an  individual  microstructure,  without  significant  changes  to 
the  basic  model.  In  addition,  a  computational  model  can  most  likely  be  developed  on  common  commercial 
software  more  quickly  than  analytic  problems  of  this  complexity  can  be  solved. 

Measuring  temperatures  experimentally  may  also  be  useful  in  predicting  problems  with  the  exposure 
process.  This  method  limits  data  collection  to  a  small  number  of  points,  and  is  not  very  useful  for 
calculation  of  deformations.  However,  the  experiment  described  in  Chapter  5  was  simple  and  inexpensive, 
and  took  little  time  to  perform.  Careful  experimentation  may  be  the  best  method  of  determining  the 
temperature  ranges  to  be  expected  during  exposure,  when  a  detailed  analysis  of  deformations  is  not 


required. 


APPENDIX 


COMPUTER  CODE 

All  of  the  code  used  for  evaluations  of  the  solutions  in  this  thesis  was  written  in  C++,  and  compiled  on 
Borland  C++  for  Windows.  Proper  allignment  of  the  characters  from  one  line  to  the  next  makes  the  code 

easier  to  read  and  understand,  so  a  courier  font  is  used;  this  is  the  default  font  in  the  Borland  C  and  C++ 
editors. 


Steady  State  Part  of  Single  Laver  Solution 
/*  file:  l“dss.cpp 

This  program  calculates  the  temperature  rise  in  a  two  dimensional 
slab  with  nonhomogeneous  boundary  conditions  on  the  upper  and  lower 
surfaces  (convection  at  y=l,  prescribed  temperature  at  y=0.)  The  ends 
are  at  zero  temperature  rise. 

*/ 

#include<math. h> 

#include<iostreain.  h> 

#include<iomanip .h> 

#include<stdio . h> 

#include<conio . h> 

#include<stdlib.h> 

tdefine  NMAX  1120  //  maximum  iterations  for  sum 

#define  PI  (4*atan(l.)) 

#define  XPOINTS  10.  //  points  in  data  grid 

#define  YPOINTS  10. 

#define  outputfile  '•2-dss.csV 
variables 

width  of  resist (x  direction) 
thermal  diffusivity 
eigenvalues  in  the  x  direction 
convection  coefficient  on  upper  surface 
thermal  conductivity 
thickness  of  resist 
summation  variable  variable 
summation  of  term  in  the  series 

constants  for  1-d  solution 

constants  for  2-d  solution 
h/k;  a  simplification 


a 

alpha 

eta 

h 

k 

L 

n 

sum 

A 

B 

C 

D 

H 


67 


truncation  value  of  infinite  sum  *-  x  direction 
temp  on  lower  surface 
temp  on  upper  surface 
number  of  data  points  in  x  direction 
number  of  data  points  in  y  direction 

//  define  global  constants 
void  main (void)  { 

/ /  initialize  constants 
double  L=0.01;  // 

double  k=0. 00198;  // 

double  h=0.16;  // 

double  a=6.0;  // 

double  Tl=2;  // 

double  Tu=5;  // 

double  H=h/k; 

double  s\im,x,y,xincr,yincr; 
double  A,B,C,D; 
double  eta; 
int  n; 

xincr=a/XPOINTS; 
yincr=L/ YPOINTS ; 

//  open  output  file 
FILE  ^outfile; 
outfile=fopen (outputf ile,  "wt")  ; 
fprintf  (outfile,  ’’dummy,  x — >\n'’ )  ; 

/ /  calculate  temperatures 
cout« "  Summation "  «endl  ; 
for (x=0 . 0 ;x<=a;x=x+xincr)  { 

f or (y=0 ; y<= (L+yincr/2 ) ;y=y+yincr)  { 
sum=0 . 0 ; 

for (n=l;n<=NMAX;n++)  { 
eta=n*PI/a; 

C=2* (l-pow{-l,n) ) /n/PI* (H*Tu-Tl* (eta*sinh (eta*L) +H*cosh (eta*L) ) ) / 
(eta*cosh(eta*L) +H*sinh(eta*L) ) ; 

D=2* (l-pow(-l,n) ) *Tl/n/PI; 

sum=siim+sin(eta*x)  *  (C*sinh {eta*y) +D*cosh  (eta*y)  )  ; 

}  //  end  n 

fprintf  (outfile,  ’•,%G",sum)  ; 
if  (  (x==a/2)£c&{y==0)  ) 

cout<< "bottom  :  *'<<sum<<endl ; 

}  //  end  y 
if ( (x==a/2) ) 

cout<<"top  :  "<<sum<<endl ; 
fprintf (outfile, "\n" ) ; 

}  //  end  X 
fclose (outfile) ; 


cm 

watts/cm/K 

watt/cm2/K 

cm 

K 

k 


NMAX 

T1 

Tu 

XPOINTS 

YPOINTS 

*/ 


//  calculate  correspoding  one-d  solution  at  top  and  bottom 
cout«endl«“One  -  d  solution " «endl  ; 

B=T1  ; 

A=H* (Tu-Tl) / {l.+L*H) ; 
cout« "bottom  :  ''«B«endl; 
cout«"top  :  "«(A*L+B)«endl; 

}  //  end  main 


Eigenvalue  Functions 
/*  file:  egnvalus.cpp 

This  program  calculates  the  eigenvalue  function  as  a  function  of 
lambda,  producing  data  for  a  graph.  Roots  are  the  eigenvalues. 
*/ 

#include<math . h> 

#include<iostream. h> 

# inc lude< i omanip . h> 

# inc lude<s  tdio . h> 

#include<conio . h> 

#include<stdlib.h> 

#include<fstream.h> 

#include<cstring.h> 

#define  RANGE  700.0 
tdefine  POINTS  1000 
#define  INPUTFILE  "solutn_g.dat" 

# define  OUTPUTFILE  " evgraph . csv" 

I*  variables 

elements  of  the  eigenvalue  matrix 
thermal  diffusivity 
value  of  determinant 
preliminary  calculation  result 
preliminary  calc,  result 
heat  transfer  coefficient 
thermal  conductivity 
distance  between  data  points 
h/k;  a  simplification 
thickness  of  the  resist 
thickness  of  the  substrate 
number  of  data  points 
range  of  lambda  for  calculations 
substrate  ambient  temperature 
resist  ambient  temperature 

main (void)  { 

double  all,al2,al3,al4; 
double  a21,a22,a23,a24; 
double  a31,a32,a33,a34; 
double  a41,a42,a43,a44; 

double  Ll,L2,hl,h2,h3,kl,k2,K,Hl,H2,H3,g,e; 
double  alphal,alpha2,det; 
double  a,Tu,Tl,step; 


a 

alpha 

det 

e 

g 

h 

k 

step 

H 

LI 

L2 

POINTS 

RANGE 

T1 

Tu 

*/ 


string  dummy; 

//  open  input  file 
f stream  input; 

input . open ( INPUTFILE, ios : : in) ; 
get line ( input , dummy) ; 

cout<<dummy<<endl<<  •* - "  <<endl  ; 

get line ( input , dummy) ; 
getline (input, dummy) ; 

//  read  values  from  input  file 
input»Ll; 

c out << dummy << "  "<<Ll<<endl; 

getline ( input , dummy) ; 
input»L2  ; 

cout<<dummy<< "  "<<L2<<endl; 

getline ( input , dummy) ; 
input>>a; 

cout«dummy«"  ''<<a«endl; 

getline ( input , dummy) ; 
input>>hl ; 

cout<<dummy« "  "<<hl<<endl; 

getline (input, dummy) ; 
input>>Tu; 

cout<<dummy<< "  "<<Tu<<endl; 

getline (input, dummy) ; 
input>>h2 ; 

c  out<<  dummy  << "  '‘<<h2<<endl; 

getline (input, dummy) ; 
input>>h3 ; 

cout<<dummy<< "  "<<h3<<endl; 

getline (input, dummy) ; 
input>>Tl ; 

cout<<dummy<< "  ''<<Tl<<endl; 

getline (input, dummy) ; 
input>>kl; 

cout<<duramy<<  *  ''<<kl<<endl; 

getline ( input , dummy) ; 
input »k2  ; 

c  ou  t «  dummy «"  *<<k2<<endl; 

getline (input, dummy) ; 
input>>alphal ; 

cout<<dummy<< "  '•<<alphal«endl; 

getline (input, dummy) ; 
input>>alpha2 ; 

cout<<dummy<< "  '•<<alpha2<<endl  ; 

//  close  input  file 
input. close ( ) ; 

//  preliminary  calculations 

K  =  kl/k2;  //  ratio  of  thermal  conductivities 

HI  =  hl/kl;  //  h/k  for  front  surface 

H2  =  h2/kl;  //  interfacial  resistance  term 

H3  =  h3/k2;  If  h/k  for  back  surface 

//  open  output  file 
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FILE  * out file; 

outfile=fopen(OUTPUTFILE, “wt" ) ; 

//  loop  on  Icimbda  to  get  values  of  determinant  for  graph 
s tep=RANGE/ POINTS ; 

for  (double  laiTibda=step;  lambda<= RANGE;  lambda = lambdas  step)  { 

//  calculate  elements  of  coefficient  matrix 
//  first  calculate  simplifying  terms: 
g  =  lambda/ sqrt (alphal) ; 
e  =  lambda/ sqrt (alpha2 ) ; 

//  now  get  the  elements  themselves 

all  =  g; 

al2  =  -HI; 

al3  =  0.0; 

al4  =  0.0; 

a21  =  (-g^cos (g*Ll) ) /H2-sin (g*Ll) ; 

a22  =  -cos (g*Ll) + (g^sin (g*Ll) ) /H2 ; 

a23  =  sin(e*Ll); 

a24  =  cos (e*Ll) ; 

a31  =  K*g*cos (g*Ll )  ; 

a32  =  0 . 0-K*g*sin (g^Ll)  ; 

a33  =  0 . 0-e*cos (e*Ll) ; 

a34  =  e*sin(e*Ll); 

a41  =  0.0; 

a42  =  0.0; 

a43  =  e*cos (e*L2 ) /H3+sin (e*L2 ) ; 
a44  =  -e^sin  (e’^L2 ) /H3  +  COS  (e*L2  )  ; 

//  find  the  determinant  of  the  matrix 

/*  det=all* (a22* (a33*a44-a34*a43 ) -a23^ (a32*a44-a34*a42 ) +a24* (a32*a43-a33*a42) ) 
-al2* (a21^ (a33*a44-a34*a43 ) -a23* (a31*a44-a34*a41) +a24* (a31*a43-a33 *a41 ) ) 
+al3* (a21* (a32*a44-a34*a42 ) -a22* (a31*a44-a34*a41 ) +a24* (a31*a42-a32*a41) ) 
-al4^ (a21* {a32*a43-a33*a42 ) -a22* (a31*a43-a33*a41) +a23* (a31*a42-a32*a41) ) ; 

*/ 

//  output  term  for  this  value  of  lambda 
det=Hl+lambda/ tan (lambda* LI) ; 
fprintf  (out file,  "%G'' ,  lambda)  ; 
fprintf  (outf  ile,  " ,  %G'‘ ,  det)  ; 
fprintf  (outf  ile,  •*\n*' )  ; 

}  //  end  for (lambda) 
fclose (outf ile) ; 

)  II  end  main 


Superposition  Solution 
/*  program  thesis!. cpp 

This  program  calculates  the  temperature  increase  in  a  slab  with  a  moving 
energy  source.  The  source  is  exponential,  and  the  slab  is  insulated  on  the 
left  and  right  edges  with  t  =  0  on  the  bottom  surface  cind  convection  on  the 
top  surface.  Multiple  passes  are  considered;  data  is  temp  as  a  function  of 
y  and  time . * / 

# include<math , h> 

# inc lude< i os  tr earn . h> 
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#include<f stream. h> 
#include<cs tring . h> 
#include<iomanip . h> 
#include<stdio .h> 
#include<conio . h> 

# inc lude<s  tdl ib . h> 

# i nc 1 ude <  t ime . h> 

#define  UNUSED  (-1.11E30) 
#define  PI  4*atan(l.) 

#define  OUTFILE  " thesis3 . csv" 
tdefine  INPUTFILE  " wafer. dat“ 
#define  MAXIT  60 


/*  variables 

a  width  of  resist (x  direction) 

alpha  thermal  diffusivity 

be[n]  eigenvalues  in  the  x  direction 

bn  eigenvalue  for  present  summation  term 

bnsq  bn  squared 

dec  sum  of  decay  solution  coefficients 

g  uniform  heat  generation 

h  heat  transfer  coefficient 

i,j  loop  variables 

k  thermal  conductivity 

1 [m]  eigenvalues  in  y  direction 

Im  eigenvalue  for  present  summation 

Imsq  Im  squared 

loop  time  required  for  one  loop 

lower  lower  limit  for  root-finding  algorithm  (eigenvalue  search) 

mu  linear  absorption  coefficient 

now  present  clock  time 

npass  number  of  passes  to  current  time 

pass  time  required  for  one  pass  of  source 

remain  time  remaining  for  calculations 

rs  spacing  between  roots  for  eigenfunction  in  y 

sterm  portion  of  answer  derived  from  the  source  term 

t  time  of  last  pass 

time  maximum  elapsed  time 

tterm  portion  of  answer  derived  from  the  time  term 

upper  upper  limit  for  root-finding  algorithm 

V  velocity  of  X-ray  beam 

w  width  of  the  X-ray  beam 

X,  y  coordinate  variables 

Amn  single  pass  coefficients 

Bmn  part  of  decay  solution  coefficients 

H  h/k  -  simplification  used  in  derivation 

L  thickness  of  resist 

MAXIT  maximum  #  of  iterations  for  ridder’s  method 

MMAX  truncation  value  of  infinite  sum  -  y  direction 

N  orthogonality  constant 

NMAX  truncation  value  of  infinite  sum  -  x  direction 

W  irradiance  at  surface 

XPOINTS  number  of  data  points  in  x  direction 


YPOINTS  number  of  data  points  in  y  direction 

*/ 

//  define  global  constants 
double  L=0 . 03 ;  II  cm 

double  k=0. 00198;  //  watts/cm/K 

double  h=0.010;  //  watt/cm2/K 

double  H=h/k; 

//  declare  subroutines 
double  froot(double  y) ; 

double  ridder (double  xl, double  x2, double  xacc) ; 
void  f ileerror (void) ; 


void  main (void)  { 

//  initialize  constants 

double  sum[21] , x, y , yincr , pass , p, t,dir,decl; 
double  Imsq, lm,bn,bnsq,b212,Amn; 
double  1ml , lm2 , N , Bmn , dec , time , tincr ; 


double  XPOINTS, YPOINTS, NMAX,MMAX; 
double  1[30+1]; 
double  be(100+l]; 
double  W,mu, V, w, a, tmax, alpha; 

//  open  input  file 
f stream  input; 

input . open ( INPUTFILE , { ios : : in  |  ios : : nocreate ) ) ; 
if ( ! input)  f ileerror ( ) ; 

//  open  output  file 
fstream  log; 


log . open (OUTFILE, ios : : out)  ; 


if ( !log)  fileerror 0 ; 
string  dummy; 

//  input  constants 
log<<'’2-d;  one  layer 
get line ( input , dummy) 
get line ( input , dummy) 
get line ( input , dummy) 
getline ( input , dummy) 
get line ( input , dummy) 
getline ( input , dummy) 
getline ( input , dummy) 
getline ( input, dummy) 
getline ( input , dummy) 
getline (input, dummy) 
getline ( input , dummy) 
get 1 ine ( input , dummy ) 
getline  ( input ,  dximmy) 
getline ( input , dummy) 
getline ( input, dummy) 
getline ( input, dummy) 
getline  ( input ,  dummy ) 
getline ( input , dummy) 
getline ( input , dummy) 


moving  heat  source  (mult  passes)"«endl; 

getline (input, dummy) ; 
input>>L;  log<<'*L;  *'<<L<<endl; 

getline (input, dummy) ; 
input>>a;  log«'‘a:  ’*<<a<<endl; 
input>>h;  log«’'h:  "«h<<endl; 
getline ( input , dummy) ; 
getline ( input , dummy) ; 
getline ( input , dummy) ; 
getline ( input , dummy) ; 
input>>k;  log«''k:  ''<<k<<endl; 
getline (input , dummy) ; 

input>>alpha;  log«" alpha:  ''<<alpha«endl; 

getline (input, dummy) ; 
input>>W;  log«*'W:  "<<W«endl; 

getline ( input , dummy) ; 
input>>mu;  log«’’mu:  "<<mu<<endl; 

getline ( input , dummy) ; 
input>>tmax;  log<<*'tmax:  *'«tmax«endl; 
input»v;  log«''v:  "«v«endl; 


getline  ( input ,  duiTttny) 
get line ( input , dummy) 
get line ( input , dummy) 


input>>w;  log<<“w:  *'<<w<<endl; 
input>>NMAX;  log<<’'NMAX:  " <<NMAX«endl  ; 

^ _ _ _ ^ _ , _ ,  input»MMAX;  log«"MMAX:  "<<MMAX«endl; 

getline  (input,  dummy)  ;  input»XPOINTS ;  log<<”x  position:  '’«XPOINTS<<endl; 

ge t line  (input,  dummy)  ;  input>>YPOINTS;  log<<''y  position:  "<<YPOINTS<<endl ; 

getline  (input,  dummy)  ;  input>>tincr ;  log<<''time  incr:  ''«tincr«endl; 

input .close ( )  ; 

log«endl«endl  ; 

double  rs=PI/L; 

int  m,n,npass; 

docket  start, now ; 

yincr=L/YPOINTS; 

//  start  the  timer  clock 
start=clock ( ) ; 

//  calculate  eigenvalues  in  y  direction 
cout<< *'  calculating  eigenvalues " «endl ; 
for (int  i=l ; i<=MMAX; i++)  { 
lml=(i-.999)*rs; 
lm2=(i-,001) *rs; 

1 [i] =ridder ( 1ml , lm2 , . 00001) ; 
cout<<i«'’  '’«1  [i]  <<endl  ; 


} 

//  calculate  eigenvalues  in  x  direction 
for(int  j=l; j<=NMAX; j++)  { 

be[j]=j*PI/a; 

} 

//  set  up  output  file 
cout<<‘* 

for (y=L;y>= (O-yincr/2 ) ;y=y-yincr)  { 
log«'' ,  •’«y; 

}  //  end  X 
log<<endl ; 


//  calculate  temperatures 

x=XPOINTS ; 

pass=a“W/v; 

log«"y->"  ; 

y=0; 

for ( i=l ; i<=10 ; i++)  { 

y=i/10.0*L; 
log«'',"«y;  } 
log<<endl ; 

f or ( time=80 .0; time<tmax; time=time+tincr )  { 
log<<time; 
npass=time/pass ; 
t=time~npass*pass ; 

//  for (y=L;y>= (yincr/2) ;y=y-yincr)  { 
f or (i=l; i<=20 ; i++) 
s\im[i]  =0.0; 

for  (m=l  ;ni<=MMAX;m++)  { 
lm=l [m]  ; 
lmsq=lm*lm; 


for {n=l;n<=NMAX;n++}  { 
bn=be[n] ; 
bnsq=bn*bn; 
b2 1 2 =bnsq+ Ims q ; 

N=  (L*  (lmsq+H*H)  +H*H)  /4*a/  (lmsq+H*H)  ; 

Bmn=-alpha/N*W*mu/k/  (bn*  (mu*mu+lmsq)  )  *  (mu*sin  ( lm*L)  -lm*cos  ( lin*L)  + 
exp (-mu*L) *lm) / (pow( (alpha*b212) , 2 ) +bnsq*v*v) * (alpha*b212* 

(cos (bn*a) -cos (bn* (a-w) ) -exp ( -alpha*b212* (a-w) /v) * (cos (bn*w) -1 . ) ) 
+bn*v* (sin (bn* a) -sin (bn* (a-w) ) -exp ( -alpha *b2 12  * (a-w) /v) * 
sin(bn*w) ) ) ; 
dec=:0 . 0  ; 
decl=0 . 0 ; 
dir=1.0; 

for (p=time“pass ;p>=0 ;p=p-pass) { 
if (dir<0) 

dec=dec+Bmn*exp (-alpha*b212*p) ; 
else 

decl=decl+Bmn*exp (-alpha*b212*p) ; 
dir=(-1.0)*dir; 

} 

Amn=-alpha/N*W*mu/k/  (bn*  (mu*inu+lmsq)  )  *  (itiu*sin  ( lm*L)  -lni*cos  ( lm*L)  + 
exp (-inu*L) *lm) / (pow ( (alpha*b212 )  ,2)  +bnsq*v*v) * (alpha*b212* 

(cos (bn* ( t*v+w) ) -cos (bn*t*v) -exp ( -alpha*b212*t) * (cos (bn*w) -1 . ) ) 
+bn*v* (sin (bn* ( t*v+w) ) -sin (bn* t*v) -exp ( -alpha*b212* t ) *sin (bn*w)  ) ) 
y=0.0; 

for (i=0;i<=10;i=i+l) { 

y=i/10.0*L; 
if (dir<0) 

sum[i] =sum[i] + (Amn+dec) *sin (bn*x) *sin (lm*y) +decl*sin (bn* (a- 
x) ) *sin ( lm*y) ; 

else 

sum[i)  =sum[ i]  (dec)  *sin (bn*x)  *sin  ( lm*y)  +  ( Amn+dec  1)  *sin (bn*  (a- 
x)  ) *sin ( lm*y) ; 

}  //  end  i  (y  loop) 

}  //  end  n 
}  //  end  m 
for ( i=l ; i<=10 ; i++ ) 
log<<'' ,  ''«sum[i]  ; 
log«endl; 
cout<<  time«endl  ; 

}  //  end  time 
clrscr ( ) ; 
now=clock( )  ; 

cout<<“ elapsed  time (min) :  " <<setprecision (3) «( (now-start) /CLK_TCK/60 . 0) ; 

log. close  0 ; 

}  //  end  main 

//  subroutine  to  tell  user  if  an  error  is  caught  in  the  input  file, 

//  or  is  the  file  isn't  in  the  working  directory, 
void  fileerror (void)  { 

cout<<’'File  error  check  for : '‘<<endl ; 

cout<<''  —  existence  of  wafer.dat  in  directory" «endl ; 
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cout«''  —  proper  format  of  file"; 
abort ( ) ; 

}  //  end  f ileerror 

double  froot (double  y)  { 

/ /  evaluates  root  function  for  eigenvalues 
double  value; 
value=y * 1 / tan ( y*L) +H ; 
return  value; 

}  //  end  froot 

double  ridder (double  xl, double  x2, double  xacc)  { 

//  ridder ‘s  method  from  "numerical  recipes  in  C, "  Press  et.al. 
/ /  finds  roots  of  function  in  froot  given  bracketting  values 
int  j  ; 

double  ans , fh, f 1, fm, fnew, s , xh, xl , xm, xnew; 
f l=froot (xl) ; 
fh=froot (x2 ) ; 

if  ( (fl>0.0&&fh<0.0)  I  I (fl<0.0&&fh>0.0) )  { 
xl=xl; 
xh=x2 ; 
ans=UNUSED; 

f or ( j  =1 ; j  <=MAXIT ; j  ++ )  { 
xm=0 .5* (xl+xh) ; 
fm=  froot  (xm); 
s=sqrt { fm*fm-f l*fh) ; 
if (s==0 . 0) return  ans; 

xnew=xm+ (xm-xl) * { ( f l>fh  ?  1.0:-1.0)*fm/s); 

if (fabs (xnew-ans) <=xacc)  return  ans; 

ans=xnew; 

fnew= froot (ans) ; 

if ( fnew==0 . 0 ) return  ans; 

if { ( (fnew) >0 . 0  ?  fabs ( fm) : -fabs ( fm) )! =fm)  { 
xl=xm; 
fl  =  fm; 
xh=ans ; 
fh=fnew; 

}  else  if (( (fnew) >0.0  ?fabs (fl) : -fabs (fl) ) ! =fl)  { 
xh=ans ; 
fh=fnew; 

}  else  if (( (fnew) >0.0  ?fabs ( fh) : -fabs ( fh) ) ! =fh)  { 
xl=ans ; 
f l=xnew; 

}  else  { 

cout<<" error  in  usbroutine  ridder"; 
abort ; 

} 

if (fabs (xh-xl) <=xacc)  return  ans; 

} 

cout<<"ridder  exceeded  maximum  iterations"; 
abort ( ) ; 

}  //  end  if 


else  cout«''root  must  be  bracketed  in  ridder"; 
return  0.0; 

)  //  end  ridder 


Two-Laver  Solution 


Program  thesis4.cpp 


This  program  calculates  the  temperature  increase  in  a  two  layer 
wafer.  Model  is  one-dimensional  with  nonhomogeneous  boundary 
conditions.  Output  is  1-D  temperature  profile  in  y  direction. 

-  definitions  - 

c  output  loop  counter 

f  output  loop  counter 

gamma  simplification  variables 

gterml,  gterm2  generation  terms  in  theta  part  of  solution 

i 
m 

min 

mul , mu2 
now 
nowt 
sec 
start 
stp 
t 

tm 

tstart 

V 

w 

y 


generic  loop  variable 
loop  counter 

number  of  minutes  elapsed 

absorption  coefficients  for  resist  and  substrate 
present  clock  time 

clock  ticks  elapsed  sinlce  beginning  of  program 

number  of  seconds  elapsed 

clock  time  at  start  of  calculations 

step  used  in  eigenvalue  search 

time  for  mathematical  model 

system  time 

time  at  start  of  theta  calculations 
scanning  velocity  of  wafer 
width  of  incident  radiation  beam 
spacial  variable 


Aim,  A2m 

Blm,  B2m 

Cl,  C2 

Dl,  D2 

PI 

LI 

L2 


eigenvalue  coefficients 

eigenvalue  coefficients 

coefficients  for  steady  state  solution 

coefficients  for  steady  state  solution 

duh! 

y  position  at  interface 
y  position  at  back  of  resist 


Tl/  Tu  temperatures  on  substrate  and  resist  surfaces,  respectively 

Tres[],  Tsub[]  temperatures  in  resist  and  substrate,  respectively 
Wl,  W2  irradiation  constant  -  layers  one  and  two 

*/ 

# include<math . h> 

#include<iostream.h> 

# inc lude< iomanip . h> 

#include<cstring,h> 

# include<s  tdio . h> 

#include<conio .h> 

#include<stdlib.h> 

#include  <malloc.h> 


# inc lude< f s tream . h> 

# inc lude<  t ime . h> 

#include<dos . h> 

#define  MAXIT  60 
#define  UNUSED  (-1.11E30) 

#define  INPUTFILE  “ wafer. daf 
#define  OUTPUTFILE  " two layer . csv" 

//  define  an  easy  squaring  function  for  later  use 
static  long  double  sqrarg; 

#define  SQR(num)  ( (sqrarg= (num) ) ==0 . 0  ?  0.0  :  sqrarg* sqrarg) 
double  lainm[301]  ; 

//  define  the  subroutines 
void  fileerror (void) ; 

double  eval (double  A, double  B, double  C, double  D, double  yl, double  y2 , 
double  1, double  a, double  p) ; 

double  lambda (double  LI, double  L2, double  kl, double  k2, double  alphal, 

double  alpha2 , double  K, double  HI , double  H2 , double  H3 , int  MMAX) 

void  main (void)  { 

/ / declarations 

int  i, c,m,min,sec; 

double  MMAX, RESPOINTS, SUBPOINTS, a, stp, tincr , y, bcterm, gterm; 

double  Tu , T1 , LI , L2 , hi , h2 , h3 , kl , k2 , K , HI , H2 , H3 ; 

double  alphal , alpha2 ; 

double  gamma, eta, Aim, A2m, Blm, B2m,N; 

double  W1 , W2 , mul , mu2 , t ; 

double  C1,C2,D1,D2; 

time_t  start, now; 

struct  time  tm; 

//  start  the  timer  clock 
start=time(NULL) ; 

//  open  input  file 
f stream  input; 

input . open ( INPUTFILE , ios : : in  |  ios : : nocreate ) ; 
if (! input)  fileerror; 

//  open  output  file 
f stream  log; 

log. open (OUTPUTFILE, ios: :out) ; 
if(!log)  fileerror; 
string  dummy; 

//  input  constants 

log<<''l"d;  two  layer;  data  in  y  and  time"«endl; 
get line (input, dummy) ;  ge t line (input , dummy) ; 

ge  t  line  (input ,  dummy)  ;  input>>Ll;  log«"Ll:  •«Ll«endl; 
ge  t  line  (input ,  dummy)  ;  input»L2;  log«''L2:  ’’«L2«endl; 
ge t line  (input,  dummy)  ;  input»a;  log<<''a:  ''<<a<<endl; 
ge  t  line  (input,  dummy)  ;  input»hl;  log«"hl:  "«hl«endl; 
getline ( input , dummy) ;  input>>Tu;  log<<"Tu:  "<<Tu<<endl; 
ge  t  line  (input ,  dummy)  ;  input>>h2;  log<<"h2:  •'<<h2<<endl; 
getline  (input,  dummy)  ;  input»h3;  log«"h3:  ’•«h3«endl; 
getline  (input ,  dummy)  ;  input>>Tl;  log«"Tl:  "«Tl«endl; 
getline  (input,  dummy)  ;  input»kl;  log«"kl:  "«kl«endl; 


"  «MMAX<<endl ; 


ge  t  line  (input,  dummy)  ;  input>>k2;  log<<"k2:  ’'«k2«endl; 

getline  ( input ,  dummy)  ;  input>>alphal ;  log<<  "alphal :  ''<<alphal<<endl; 

get line ( input , diommy) ;  input>>alpha2 ;  log<<"alpha2 :  "<<alpha2<<endl; 

getline  (input ,  dummy)  ;  input>>Wl;  log<<'*Wl:  ''«Wl<<endl; 

getline  (input,  dummy)  ;  input»W2;  log«"W2:  '’«W2«endl; 

getline  ( input ,  dximmy)  ;  input>>mul;  log<<’'mul:  '’<<mul<<endl; 

getline  ( input ,  d\immy)  ;  input>>mu2  ;  log«  ’*mu2  :  '•«mu2«endl  ; 

getline  (input,  dummy)  ;  input>>t;  log<<’‘tmax:  '’<<t«endl; 

getline  (input ,  d\immy)  ;  getline  ( input ,  dummy)  ; 

getline (input , dummy) ;  getline (input , dummy) ; 

getline (input, dummy) ;  getline (input , dummy) ; 

getline  ( input ,  dummy)  ;  input>>MMAX;  log<<''^u.^  " <<MMAX<<endl ; 

getline  (input,  dummy)  ;  input»RES POINTS;  log« " RESPOINTS :  ■‘«RESPOINTS«endl; 
getline  ( input ,  dummy)  ;  input>>SUBPOINTS ;  log<< "  SUBPOINTS :  **  <<SUBPOINTS<<endl ; 

getline  (input,  dummy)  ;  input»tincr ;  log«’'time  increment:  "  «tincr<<endl  ; 
log<<endl«endl ; 

//  close  input  file 
input . close ( ) ; 
double  *Tres,*Tsub; 

if  (RESPOINTS>40  1|  SUBPOINTS>40 )  { 

cout<<" Exceeded  40  data  points 
abort ( ) ; 

} 

try  { 

Tres=new  double  [41]  ; 

Tsub=new  double  [41]  ; 


catch (xalloc)  { 

cout«endl«" could  not  allocate  memory"; 
exit (-1) ; 

) 

//  initialize  tres  and  tsub 
for ( i=0 ; i<=RESPOINTS ; i++ )  { 

Tres [i] =0 . 0 ; 

) 

for (i=0; i<=SUBPOINTS; i++)  { 

Tsub [i] =0 . 0 ; 


//  preliminary  calculations 

K  =  kl/k2;  //  ratio  of  thermal  conductivities 

HI  =  hl/kl;  //  h/k  for  front  surface 

H2  =  h2/kl;  //  interfacial  resistance  term 

H3  =  h3/k2;  //  h/k  for  back  surface 

//  find  the  eigenvalues  in  the  y  direction  for  the  theta  (time 
//  dependent)  part  of  the  solution 

s  tp= lambda ( LI , L2 , kl , k2 , alphal , alpha2 , K , HI , H2 , H3 , MMAX) ; 
lamm[0] =0 ; 

cout«endl«" calculating  temperatures.  . 

//  calculate  temperatures 
f  or  (m=l  ;m<=MMAX;m++ )  { 

gammas lamm[m]  /sqrt  (alphal)  ; 
eta=lamm[m]  /  sqrt  (alpha2 )  ; 


80 


//  define  the  constants  for  the  y  direction  eigenfunctions 
A2ra=l . 0 ; 

B2m=  (-  (eta'^cosl  (eta*L2 )  /H3  +  sinl  (eta*L2 )  )  )  / 

(-eta’^sinl  (eta*L2 )  /H3+cosl  (eta*L2 )  )  ; 

Blm= (-B2m*cosl (eta*Ll) ’Sinl (eta*Ll) ) / (Hl/gamma* (-gamma*cosl (gamma* LI) /H2- 
sinl (gamma*Ll ) ) +gamma*sinl (gamma*Ll ) /H2-cosl (gamma*Ll) ) ; 

Alm=Blm*Hl  /  gamma  ; 

//  define  orthogonality  constant 

N=kl/ (alphal* gamma) * (Alm*Alra* (gamma*Ll/2 . -sinl (2 . *gamma*Ll) * .25) - 

Alm*Blm*cos (2 . *gamma*Ll) /2 .+Blm*Blm* (gamma*Ll/2 . +sin ( 2 . *gamma*Ll) /4 . ) 
+Alm*Blm/2. )  + 

k2/alpha2/eta* (A2m*A2m* (eta*L2/2 . -sin (2 . *eta*L2 ) /4 . ) 

-A2m*B2m*cos (2 . *eta*L2 ) /2 . +B2m*B2m* (eta*L2 /2 . +sin ( 2 . *eta*L2 ) /4 . ) - 
A2m*A2m* (eta*Ll/2 . -sin (2 . *eta*Ll) /4. ) +A2m*B2m*cos ( 2 . *eta*Ll ) /2 
B2m*B2m* (eta*Ll/2 . +sin (2 . *eta*Ll) /4 . )  ); 

gterm=Wl*mul/ {SQR(mul) +SQR(gamma) ) * (  Aim* (  mul*sin(gamma*Ll)- 

gamma*cos (gamma*Ll) +gamma*exp (-mul*Ll)  ) +Blm* (  mul*cos (gamma*Ll ) + 
gamma*sin(gamma*Ll) -gamma*exp (-mul*Ll)  )  ) +W2*mu2/ (SQR(mu2 ) +SQR (eta) ) * 

(  A2m* (  exp(-mu2*Ll)* (mu2*sin(eta*L2)-eta*cos (eta*L2) )“exp(-mu2*L2) * 
(mu2 *sin (eta*Ll ) -eta* cos (eta* LI) ) ) +B2m* (exp (mu2*Ll) * (mu2*cos (eta*L2 ) + 
eta*sin(eta*L2) ) -exp (-mu2*L2) * (mu2*cos (eta*Ll) -eta*sin(eta*Ll) )  )  ) ; 

bcterm= (Tu*hl*Blm+Tl*h3* (A2m*sin (eta*L2 ) +B2m*cos ( eta*L2 ) ) ) ; 

//  add  theta  term  to  resist  solution 
for (c=0;c<=RESPOINTS;c++) { 
y=c*Ll /RESPOINTS ; 

Tres [c] =Tres [c] +1 . 0/N* (Aim* sin (gamma* y) +Blm*cos (gamma*y) ) 

*  (l-exp(-SQR(lamm[m]  )  *t)  )  /  (SQR(lamm[m]  )  )  *  (bcterm+gterm)  ; 

} 

//  add  theta  term  to  substrate  solution 
for ( c=0 ; c<=SUBPOINTS ; C++ )  { 

y=C* (L2-L1) /SUBPOINTS+Ll; 

Tsub [c] =Tsub [c] +1 . 0/N* (A2m*sinl (eta*y) +B2m*cosl (eta*y) ) 

*  (0.0*expl  (  -  (SQR(lainm[m]  )  )  *t)  +bcterm)  ; 

} 

}  //  end  m  loop 
gettime  (Sctm)  ; 

if  ( tm. ti_hour>12) tm. ti_hour=tm. ti_hour-12 ; 
printfCXn  finished  at:  %2d:  %02d:  %02d'‘ , 

tm.  ti_hour ,  tm.  ti_min,  tm.  ti_sec)  ; 

//  output  the  temperatures 
//  print  column  headings  (y  values) 
log«''  " ; 

for (c=0;c<=RESPOINTS;c++) { 
y=c  *  LI /RESPOINTS ; 
log<<’' ,  ''<<y; 

} 

f or ( c=0 ; c<=SUBPOINTS ; C++ )  { 

y=c* (L2-L1) /SUBPOINTS+Ll; 
log«" ,  "«y; 

} 

log<<endl  ; 

//  output  the  temps 


for ( c=0 ; c<=RESPOINTS ; C++ ) { 
log<< " , " <<Tres [ c ] ; 

} 

for {c=0;c<=SUBPOINTS;c++)  { 
log«  “ ,  "  «Tsub  [  c  ]  ; 

} 

log«endl<< "steady  state : " «endl ; 

//  calculate  and  output  steady  state, 

//  no  generation  temperatures 
log«" 

for (c=0;c<=RESPOINTS;c++)  { 
y=c*Ll/RESPOINTS; 

Cl= (Tl-Tu) / { 1/H1+ (1/H2+L1-K*L1+K*L2+K/H3 ) ) ; 

D1=TU+C1/H1; 
log«" ,  "«(Cl*y+Dl)  ; 

} 

for ( c=0 ; c<=SUBPOINTS ; C++ )  { 
y=:c*  (L2-L1)  /SUBPOINTS+Ll; 

Cl=- (Tl-Tu) / (1/Hl- (1/H2+L1-K*L1+K*L2+K/H3) ) ; 

C2=K*C1; 

D2=T1-C1* (K*L2+K/H3) ; 
log«" ,  "«  (C2*y+D2)  ; 

} 

log<<endl; 

//  output  iegenvalues 
for {i=l;i<=MMAX; i++)  { 

log«endl<<i« " ,  "«lamm[i]  ; 

} 

log<<endl ; 

//  delete  arrays  to  clear  memory 
delete []  Tres;  delete []  Tsub; 

//  close  output  file 
log. close ( ) ; 

cout<<endl«endl«"data  written  to  file  -  processing  completed"; 
now=  t ime ( NULL ) ; 
sec=difftime (now, start) ; 
if  (sec>60)  { 
min=sec/60 ; 
sec=sec-min*  60 ; 

cout<<endl«"  elapsed  time:  "<<min<<"m  "«sec<<"s"; 

) 

else  { 

cout«endl«"  elapsed  time:  "<<sec«"s"; 

} 

}  //  end  main  program 

//  subroutine  to  tell  user  if  an  error  is  caught  in  the  input  file, 
//  or  is  the  file  isn’t  in  the  working  directory, 
void  f ileerror (void)  { 

cout<<"File  error  check  for : "<<endl; 

cout«"  —  existence  of  "«INPUTFILE«"  in  directory"«endl ; 
cout«"  —  proper  format  of  file"; 
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abort ( )  ; 

}  //  end  fileerror 

//  evaluates  a  portion  of  the  theta  part  of  the  solution 
double  evaKdouble  A, double  B, double  C, double  D, double  yl, double  y2, 
double  1, double  a, double  p)  { 
double  ret,  soa, ch2 , sh2 , c2 , s2 , chi , shl , cl, si; 
soa  =  sqrt(a); 
ch2=cosh  (p*y2)  ,- 
sh2=sinh(p*y2) ; 
c2=cos {l*y2/soa) ; 
s2=:sin{l*y2/soa)  ; 
chl=cosh (p*yl) ; 
shl=sinh (p*yl) ; 
cl=cos (l*yl/soa) ; 
sl=sin(l*yl/soa) ; 
ret=l . 0/ (p*p+l*l/a) * 

(A*C* {p*ch2*s2-l/soa*sh2*c2) +B*C* (p*ch2 *c2+l/soa*sh2*s2 ) + 

A*D* (p*sh2*s2“l/soa*ch2*c2) +B*D* {p*sh2*c2+l/soa*ch2*s2) - 
(A^C* (p*chl*sl-l/soa*shl*cl)+B*C* (p*chl^cl+l/soa*shl*sl ) + 

A*D* (p*shl*sl-l/soa*chl*cl)+B*D* (p*shl*cl+l/soa*chl*sl) ) ) ; 
return (ret) ; 

)  //  end  eval 

//  subroutine  to  calculate  the  value  of  the  deteriainant  of  the 
//  coefficient  matrix  given  a  guess  at  lambda  (roots  of  the 

//  matrix  are  the  eigenvalues  in  the  y  direction,) 

long  double  f root (double  lambda, double  alphal , double  alpha2 , double  LI, 
double  L2, double  K, double  HI, double  H2, double  H3)  { 

//  This  subroutine  finds  the  value  of  the  determinant  for  a 
if  given  value  of  lambda. 

long  double  all , al2 , al3 , al4 ; 
long  double  a21 , a22 , a23 , a24 ; 
long  double  a31 , a32 , a33 , a34 ; 
long  double  a41 , a42 , a43 , a44 ; 
long  double  g,e,det; 

/ /  calculate  elements  of  coefficient  matrix 
//  first  calculate  simplifying  terms: 
g  =  lambda/ sqrt (alphal) ; 
e  =  lambda/ sqrt (alpha2) ; 

/ /  now  get  the  elements  themselves 

all  =  g; 

al2  =  -HI; 

al3  =  0.0; 

al4  =  0.0; 

a21  =  g*cos(g*Ll) /H2+sin(g*Ll) ; 

a22  =  cos(g*Ll)-g*sin(g*Ll)/H2; 

a23  =  ~sin(e*Ll); 

a24  =  -cos(e*Ll); 

a31  =  K*g*cos (g*Ll) ; 

a32  =  0. 0-K*g*sin(g*Ll) ; 

a33  =  0 . 0-e*cos (e*Ll ) ; 
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a34  =  e*sin(e*Ll); 
a41  =  0.0; 
a42  =  0.0; 

a43  =  e*cos (e*L2 ) /H3+sin (e*L2 ) ; 
a44  =  -e*sin(e*L2) /H3+COS (e*L2) ; 

//  find  the  determinant  of  the  matrix 

det  =  all^ (a22* {a33*a44-a34*a43 ) -a23* (a32*a44-‘a34*a42 ) +a24* (a32*a43-a33*a42 ) ) 
-  al2* (a21* (a33*a44-a34*a43 ) -a23* (a31^a44-a34*a41) +a24* (a31*a43-a33 *a41 ) ) 

•f  al3*  (a21*  (a32^a44-a34*a42 )  -a22*  (a31*a44-a34*a41)  4.a24*  (a31*a42-a32*a41 )  ) 
al4* (a2X* (a32*a43-a33^a42 ) -a22* (a31*a43-a33*a41) +a23* (a31*a42-a32 *a41 ) ) ; 
//  return  value  of  determinant 
return  det; 

}  //  end  subroutine  froot 

//  find  roots  given  bracketing  values 

double  ridder (double  xl, double  x2, double  xacc, double  alphal , double  alpha2 , 
double  LI, double  L2 , double  K, double  HI, double  H2, double  H3)  { 

//  ridder 's  method  from  "numerical  recipes  in  C, "  Press  et.al. 

/ /  finds  roots  of  function  in  froot  given  bracketting  values 
//  xacc  is  the  accuracy  of  the  routine 
int  j  ; 

long  double  ans , fh, f 1 , fm, fnew, s , xh, xl , xm, xnew; 
f l=f root (xl , alphal , alpha2 , LI , L2 , K, HI , H2 , H3 ) ; 
f h=f root ( x2 , alphal , alpha2 , LI , L2 , K , HI , H2 , H3 ) ; 
if  ( (fl>0.0&&fh<0.0)  I  I  (fl<0.0&&fh>0.0) )  { 
xl=xl; 
xh=x2 ; 
ans=UNUSED; 

for ( j=l; j<=MAXIT; j++)  { 

xm=0 . 5* (xl+xh) ; 

fm=froot (xm, alphal , alpha2 , LI , L2 , K, HI, H2 , H3 ) ; 
s=sqrt  ( fm*fm-f l*fh)  ; 
if (s==0 . 0) return  ans; 

xnew=xm+ (xm-xl) * ( ( fl>fh  ?  1.0:-1.0)*fm/s); 
if (fabs (xnew-ans) <=xacc)  return  ans; 
ans=xnew; 

fnew=f root (ans , alphal , alpha2 , LI , L2 , K, HI , H2 , H3 ) ; 
if ( fnew==0 . 0 ) return  ans; 

if (( (fnew) >0.0  ?  fabs ( fm) ; -fabs ( fm) )! =fm)  { 
xl=xm; 
fl=fm; 
xh=:ans  ; 
fh=fnew; 

}  else  if (( (fnew) >0.0  ?fabs ( f 1) : -fabs ( f 1) ) ! =f 1)  { 
xh=ans ; 
fh=fnew; 

}  else  if (( (fnew) >0.0  ? fabs ( fh) : -fabs ( fh) )! =fh)  { 
xl=ans ; 
f l=xnew; 

}  else  { 

cout<<" error  in  subroutine  ridder"; 
abort ; 


} 

if ( fabs (xh~xl) <=xacc)  return  ans; 


} 

cout«"ridder  exceeded  maximum  iterations"; 
abort ( ) ; 

} 

else  cout«"root  must  be  bracketed  in  ridder"; 
return  0.0; 

}  //  end  ridder 

//  main  eigenvalue  finding  routine 

double  lambda (double  LI, double  L2, double  kl, double  k2, double  alphal, 
double  alpha2 , double  K, double  HI, double  H2, double  H3 , int  MMAX) 

{ 

double  rl,step,v; 
int  m,chk,inc; 
double  acc; 

/*  find  sign  changes  in  the  value  of  the  determinant  by  walking 
through  potential  values  of  lambda,  and  use  ridder 's  method  on 
the  intervals  with  sign  changes  to  find  the  roots  (eigenvalues) 

_ variable  list _ 

acc  desired  accuracy  of  the  eigenvalues 

rl  start  value  for  root  search  interval 

r2  end  value  for  interval 

m  subscript  for  lambda 

step  size  of  intervals  to  be  searched 

passed  variables  are  the  same  as  their  counterparts 
in  the  main  program. 

*/ 

step=  0.00001; 
acc  =  0.0000001; 
rl  =  0; 

//  find  first  root  with  very  small  steps 
clrscr ( ) ; 

cout<<endl« "Looking  for  first  eigenvalue.  .  ."<<endl; 
do  { 

rl-rl+step; 

} 

while ( f root ( r 1 , alphal , alpha2 , LI , L2 , K , HI , H2 , H3 ) / 

froot( (rl+step) , alphal, alpha2, LI, L2,K, HI, H2,H3) >0.0) ; 

//  first  root  is  between  rl  and  rl+step;  now  find  it 

lamm [ 1 ] =r idder ( r 1 ,  ( r 1 +s  tep ) , acc , alphal , alpha2 , LI , L2 , K , HI , H2 , H3 ) ; 

step=lamm[l] /lOO . ; 

rl=lamm[l] +step; 

m=2  ; 

//  now  find  rest  of  roots  with  new  step  size 

cout<< "  searching  for  eigenvalue  2  of  "«MMAX; 

chk=0; 

inc=0; 

do  { 

v=  f root ( r 1 , alphal , alpha2 , LI , L2 , K , HI , H2 , H3 ) / 


froot ( (rl+step) ,alphal, alpha2,Ll,L2 ,K,H1,H2,H3) ; 
if(v<0.0)  { 

//  root  is  in  range;  find  it 

lamm [m] =ridder ( rl , ( rl+step ) , acc , alphal , alpha2 , LI , L2 , K , HI , H2 , H3 ) ; 

//  reset  starting  value  to  just  past  the  last  root 

rl=laiTuu[m]  +step/10 .  ; 

in++ ; 

chk=0 ; 

clrscr ( ) ; 

if  (m<=MMAX) 

cout<<'’ searching  for  eigenvalue  *'<<m<<'*  of  '‘<<MMAX<<"  above  lam="<<rl 
<<endl«'' stepping  at  ''«step<<";  increased  step  ''<<inc«"  times"; 
}  //  end  if 

else  { 

//  root  is  not  in  range 
if  (chk>3000)  { 
inc=inc+l ; 
step=step*1.4; 
chk=0; 

} 

chk=chk+l ; 
rl=rl+step; 

}  //  end  else 

}  //  end  do 
while  (m<=MMAX)  ; 
return (step) ; 

}  //  end  subroutine  lambda 


Input  File  for  Programs 

Data  for  wafer  heating  problems 


resist  thickness: 

0.03 

wafer  thickness: 

0.1 

wafer  length: 

100.0 

hi: 

0.001 

resist  surface  temp: 
10.0 
h2: 

1. 0e7 
h3: 

1.0e7 

substrate  surface  temp: 

0.0 

res.  conductivity: 
0.00198 

sub.  conductivity: 


1000 

res.  diffusivity: 

0.00001182 

sub.  diffusivity: 

1.0 

Wl: 

-10.82 
W2  : 

0.0 

mul : 

-23.10 
mu2  : 

-1.0 
time : 

4.0 

beam  velocity: 

0.01 

beam  width: 

95.0 

NMAX  (<100) 

70 

MMAX  (<30,  or  300  for  two  layers) 
29 

xpointS/  X  value,  respoints 
4 

ypoints,  subpoints 
10 

time  increment 
0.5 
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